1 Introduction

Frozen soil (which comprises permafrost and seasonally frozen soil) plays an important role in the land surface hydrology of cold regions. Freeze-thaw processes significantly modulate the hydraulic and thermal characteristics of soil. Frost and thaw fronts (FTFs) separate a soil column into frozen and unfrozen segments, and the changes in the FTFs will affect the water and energy cycles between the land surface and the atmosphere.

Permafrost (the ground that remains at or below 0 ℃ for at least two consecutive years) and seasonally frozen ground areas account for more than 30% of the Earth’s total land area^[1]. In particular, the area with frozen soil is equivalent to 75% of the land surface in China, and the permafrost is distributed mainly on the Qinghai-Tibet Plateau and in the northeast of China. The maximum frost depth is an important parameter when roadbeds in seasonally frozen areas are designed^[2-5]. Frozen soils have a limited capacity for infiltration which leads to high levels of runoff during the spring snowmelt^[6-7]. Moreover, the thaw depth of the soil is important for cultivation and seeding during early spring^[8]. Changes in the depths of the FTFs can affect ecosystems via carbon cycling. Thus, the organic matter decomposition rate in the unfrozen soil is ten times larger than that in frozen soils in the black spruce forest^[9]. The increased respiration rate of soil can be explained by the increases in the thaw depth during mid-summer. The FTF depth is sensitive to the climate change, and it is an important indicator of the climate change^[10].

Interpolating the soil temperature (observed or simulated) and obtaining the position of the freezing and thawing point Tf isotherm are a typical method for estimating the depth of FTFs. However, multiple FTFs cannot be simulated in the same soil layer by use of this method. In addition, this method will produce numerical oscillations when the soil temperature is close to the freezing temperature Tf during the autumn freezing and spring snowmelt periods.

The Stefan equation^[11] can be used to calculate the FTFs in a one-directional soil column, where this method assumes that all of the heat reaching the FTFs can be used for freezing the soil water or for melting the soil ice. Little site-specific information is available. Therefore, the Stefan solution can be used frequently by permafrost scientists to predict the frost depth and to simulate the heat transfer during water phase transitions in the frozen soil^[12-14]. Woo et al.^[15] used the surface and deep soil temperature as forcing variables in the Stefan equation, and proposed a two-directional freeze and thaw algorithm.

Significant studies of soil freezing and thawing processes have been performed since the 1970s when Harlan^[16] constructed the first hydrothermal coupled-migration model. Based on field observations and laboratory experiments, soil water and heat coupling models as well as their corresponding numerical simulation algorithms have been developed for different needs^[17-21]. These studies have advanced the study of soil hydrothermal processes. However, they have not considered changes in the depth of FTFs and their influence on the model.

In this study, we use a two-directional freeze and thaw algorithm to simulate changes in the FTF depth, where this method is incorporated into the thermal diffusion equation. In the version 4.5 of the Community Land Model (CLM), the thermal diffusion equation is used to simulate the temperature for a 15-layer soil column. A local adaptive variable-grid method is used to discretize the model. The results of sensitivity tests demonstrate that this method is stable and the FTFs can be tracked continuously. The present simulation of the QinghaiTibet Plateau D66 site shows that the incorporated model performs much better in the soil temperature simulation than the original thermal diffusion equation. The soil temperature profile within a soil layer can be estimated according to the position of the freezing front, thereby demonstrating the potential application of this method to land-surface process modeling.

2 Thermal diffusion model including changes of FTFs 2.1 Thermal diffusion equation

We consider the simulation of a one-dimensional vertical soil column. Let (0, L) be the one-dimensional vertical soil column, z = 0 be the top of the soil column, and z = L be the bottom of the soil column. We assume that the position of the freezing and thawing point T_f isotherm is defined as the position of the soil FTF.

Based on the energy conservation, the thermal diffusion equation can be written as follows:

(1)

(2)

where t (s) is the time, z (m) is the space coordinate (positive downward), C (J·m⁻³·K⁻¹) is the soil heat capacity, λ (W·m⁻¹·K⁻¹) is the soil thermal conductivity, T (K) represents the soil temperature, and F (W·m⁻²) is the heat flux.

The soil heat capacity C (J·m⁻³·K⁻¹) is calculated as follows:

(3)

where θ_s is the soil porosity, and C_g, C_w, and C_i represent the volumetric heat capacities of soil solid, water, and ice, respectively.

In Eq. (3), C_g can be calculated as

(4)

where P_s and P_c are the percentages of the sand and the clay in the soil, respectively.

The soil thermal conductivity λ can be calculated as

(5)

where Ke is the Kersten number, which can be calculated as follows:

(6)

where

(7)

The saturated soil thermal conductivity λsat and the dry soil thermal conductivity λdry can be calculated as

(8)

(9)

where λ_s, λ_i, and λ_l are the thermal conductivities of soil solid, ice, and water, respectively, ρ_d (kg·m⁻³) is the bulk density, and Q is the volumetric fraction of the soil moisture content.

In Eq. (8), the thermal conductivity of soil solid λs can be calculated as

(10)

We use the following equation to calculate the supercooled water content when the soil temperature is less than T_f:

(11)

where θ_max is the maximum liquid water when the soil temperature T is below the freezing temperature T_f, g is the acceleration due to gravity, ψ_s is the saturated water potential, L_i is the latent heat of fusion, and b is the Hornberger exponent.

We assume that the heat flux is zero for the lower boundary conditions. The discretization of Eqs. (1) and (2) can be calculated as

(12)

(13)

where the parameter α = 0.5.

2.2 Calculation of FTFs

We use a two-directional freeze and thaw algorithm^{[15, 22]} to simulate the depth of the FTFs in the present study. First, we consider the case of downward freezing or melting.

The Stefan equation can be described in the following form:

(14)

(15)

where z_l (l = f or t) (m) is the frost/thaw front depth, D (℃·s) is the freeze/thaw index, t is the freezing or thawing duration, T (℃) is the average temperature at the surface, and L (J·m⁻³) is the volumetric latent heat of fusion.

Squaring the Stefan equation (14), we have

(16)

and differentiating Eq. (15) for D, we have

(17)

(18)

We assume that the thickness of the ith soil layer is Δz_i, the depth of the ith soil layer is z_i, the requirement of the freeze/thaw index for the frost/thaw front to pass from z_i−1 to z_i is N_i, θ_i is the volumetric fraction of the ith layer soil moisture, and the thermal resistance of the ith layer soil is R(i) = Δz_i/λ_i.

If z_f = 0 when t = 0, then for the first soil layer, ΔD = N₁, Δz_l = Δz₁, and z_l changes from 0 to z₁. If we divide the first soil layer into N layers, then N₁ can be calculated as follows:

(19)

When N → ∞,

(20)

analogous to this rule,

(21)

When the freeze/thaw index and ,

(22)

where z_f0 is the FTF depth in the ith soil layer.

By solving Eq. (21), we can obtain z_f0 as

(23)

and the FTF depth is

(24)

If z_l > 0, and z_i−1 < z_l ≤ z_i, when t = 0, then the freeze/thaw index can be used as follows to calculate the FTF depth:

(25)

The same approach is used for the upward calculations.

2.3 Numerical algorithm incorporating FTFs into thermal diffusion equation

If we assume that the soil temperature and the FTF depth are given at the time tk, then the values of these variables at tk+1 can be calculated as follows^[23-25]:

(i) Divide the soil column into n layers. The depth of the ith layer soil is z_i (i = 1, 2,···, n). The thickness of the ith layer soil is Δz_i (i = 1, 2,···, n). The soil heat capacity is C_i^k (i = 1 , 2,···, n), and the soil thermal conductivity is λ_i^k (i = 1, 2,···, n). Assume that the temperature of the soil layers at the time t_k is T_i^k (i = 1, 2,···, n). The number of the frost depths is J₁, and the frost depth value is z_{f, j1}^k (j₁ = 1, 2,···, J₁). The number of the thaw depths is J₂ , and the thaw depth value is z_{t, j2}^k (j₂ = 1, 2,···, J₂). If we assume that N FTFs satisfy the following conditions: z_{f, j1}^k − z_i > 0.02Δz_i and z_i+1 − z_{f, j1}^k > 0.02Δz_i, or z_{t, j2}^k − z_i > 0.02Δz_i and z_i+1 − z_{t, j2}^k > 0.02Δz_i, then we can obtain the initial values of the FTF depths by interpolating the initial soil temperature.

(ii) Add the number of nodes from n to n + N. Update z_i, Δz_i, and T_i^k to z_m, Δz_m, and T_m^k (m = 1, 2,···, n + N), respectively. If z_m = z_{f, j1}^k or z_m = z_{t, j2}^k, then T_m^k = 0 ℃. Interpolate the soil heat capacity C_i^k and the soil thermal conductivity λ_i^k to obtain C_m^k and λ_m^k.

(iii) As described in Subsection 2.1, calculate the soil temperature T_m^k+1 at the time t_k+1 for n + N nodes with Eqs. (12) and (13). According to the corresponding relationship for z_i (i = 1, 2,···, n) and z_m (m = 1, 2,···, n + N), return the updated variable soil temperature which is from T_m^k+1 (m = 1, 2,···, n + N) to T_i^k+1 (i = 1, 2,···, n).

(iv) According to Subsection 2.2, calculate the frost front depth z_{f, j1}^k+1 (j₁ = 1, 2,···, J₁) and the thaw front depth z_{t, j2}^k+1 (j₂ = 1, 2,···, J₂) at the time t_k+1 with Eqs. (22), (23), and (24). If |z_{f, j1}^k+1 − z_{t, j2}^k+1| ≤ 0.01, then z_{f, j1}^k+1 = 0, and z_{t, j2}^k+1 = 0.

(v) Return to repeat the step (i) to the step (iv).

3 Numerical experiments 3.1 Ideal test

To test and verify the correctness and numerical stability of the algorithm, the following sensitivity tests are performed. First, we compare the results of FTF simulations using two hypothetical soil profiles. Second, the FTF depths are simulated with different time steps and spatial resolutions. The details of the tests are described in the following.

3.1.1 Simulation of FTFs in multi-layered soil

The physical parameters of the two soil profiles are given in Table 1. Both of the profiles contain two soil layers, where the depths of the layers are 0.6 m and 1.5 m. Figure 1(a) is a homogeneous silt soil, and its properties remain constant with the depth. In Fig. 1(b), the upper layer is the silt soil, and the lower soil is the peat soil, which has a high water content and low thermal conductivity. In both cases, the surface temperature is the same, i.e., T_surf = −10 ℃. The initial frost depth values are both z_f = 0, and the time step is Δt = 3 600 s. We only consider the case of downward ground freezing.

Figure 1 shows the simulation results for the frost front depth in the two soil profiles. In both cases, the initial results are as expected. However, as the frost front moves through the interface between the two different layers, the frost front depth is deeper in Case 1. There is an abrupt change in the gradient in Case 2, which does not occur in Case 1. This is because the peat soil in Case 2 has a higher water content. Therefore, a lower temperature is required in the peat soil layer.

In the common form of the Stefan equation, the thermal conductivity λ and the volumetric fraction of soil moisture θ are constants. Thus, this method is only suitable for the homogeneous soil. This test shows that the new algorithm reflects the physical properties of different soils. Based on these simulations, we conclude that our algorithm can be used for multi-layered soils.

3.1.2 Comparative test

We assume that the upper boundary condition is given by the following periodic temperature boundary condition:

(1)

The initial values of the FTFs are z_f = z_t = 0. The soil hydraulic parameters are λ_l = 0.57 W·m⁻¹·K⁻¹, λ_i = 2.29 W·m⁻¹·K⁻¹, and λ_s = 7.935 3 W·m⁻¹·K^−1[26]. The simulated volumetric water content is 0.3. The soil column, which is based on the hierarchical structure of the CLM^[27], is discretized into 15 layers, where the depths of the nodes are

(2)

and f_s = 0.025 is a scaling factor.

Under the experimental conditions described above, the following control tests are performed:

Test 1 The time steps are 0.5 h, 1 h, and 2 h, respectively. The spatial resolution is based on the hierarchical structure of the CLM. The FTFs are obtained as described in Subsection 2.2.

Test 2 The spatial resolution is based on the hierarchical structure of the CLM. The spatial resolution is set as a thickness of 1 cm for each soil layer. The time step is 1 h. We then simulate the FTF depths as described in Subsection 2.2.

Figure 2 shows the calculated results for different time steps and different spatial resolutions. Figures 2(a), 2(b), and 2(c) show the simulated results for the FTF depths at the three time steps. Using the results obtained at 1 h as a reference, the error in the calculated results at the other two time steps is shown in Fig. 2(d). The maximum error in the results at 0.5 h is 0.008 m, and the maximum error in the results at 2 h is 0.018 m. Figure 2(e) shows the simulated FTF depths at the high spatial resolution (the thickness of each soil layer is 1 cm). According to Fig. 2(f), the error in the results at different spatial resolutions is less than 0.006. These results demonstrate that the numerical algorithm is correct and stable.

3.2 Simulation test for Tibetan Plateau D66 site

The D66 site (35°31′N, 93°47′E) is located in the northern part of the Tibetan Plateau at an altitude of 4 560 m, where the precipitation is low, and the type of frozen soil is permafrost. The soil comprises 58% sand and 10% clay^{[19, 28-29]}. The observations obtained from the automatic weather station at D66 include atmospheric forcing fields at 1.5 m every 30 min. The soil temperatures in 10 layers (4 cm, 20 cm, 40 cm, 60 cm, 80 cm, 100 cm, 130 cm, 160 cm, 200 cm, and 250 cm) are obtained once each hour.

To verify the superior performance of the thermal diffusion model including the changes in the FTF depth proposed in this study, the following two models are tested:

Model 1 As described in Subsection 2.3, we incorporate the two-directional freeze and thaw algorithm into the thermal diffusion equation.

Model 2 As described in Subsection 2.1, we simulate the soil temperature using the original thermal diffusion equation.

In order to facilitate their comparison, the time step, the spatial resolution, the initial condition, and the boundary condition are all the same. The depth of the soil column is 2.5 m. The node depths are 4 cm, 20 cm, 40 cm, 60 cm, 80 cm, 100 cm, 130 cm, 160 cm, 200 cm, and 250 cm, and the time step is 1 h. The upper boundary condition and the initial soil temperature value are determined based on observations, and the lower boundary condition is assumed for the zero heat flux. The simulated time period ranges from September 1, 1997 to September 22, 1998.

Figure 3 shows the FTF depth changes during the simulation period according to Subsection 2.2, which we also compare with the observations obtained at the same time. The frost duration in the surface soil is up to about 5 months. As shown in Fig. 3, the freezing process occurs over about 3 months, and it reaches 250 cm in late December. The maximum observed depth is 250 cm, therefore, the maximum FTF depth is 250 cm in Fig. 3. The soil begins to thaw in early April, and the continuous thaw reaches down to 250 cm after almost 3 months. The correlation coefficients (CCs) for the frost front and thaw front in the simulations and observations are 0.79 and 0.80, respectively, and the root mean squared errors (RMSEs) are 0.58 and 0.42.

Figure 4 shows the soil temperatures at 20 cm, 60 cm, 100 cm, and 160 cm according to the simulations of Model 1 and Model 2. Comparisons with the observations are also performed for the same time period. According to Fig. 4, the trends in the soil temperature simulated by the two models are consistent with the observations. The red lines representing Model 1 are closer to the actual values than the blue lines representing Model 2. Table 2 shows the CCs and RMSEs for the soil temperature in different soil layers using simulations of Model 1 and Model 2 and based on the observations. According to Table 2, Model 1 has lower RMSEs at depths of 20 cm, 60 cm, 100 cm, and 160 cm. Figure 4 and Table 2 demonstrate that the simulation results obtained by Model 1 are better than those using Model 2. The improvement using Model 1 compared with Model 2 is better in the lower levels than that in the upper level. This is because when the model considers the change in the FTFs, the soil temperature profile within a soil layer can be estimated according to the positions of the FTFs. Comprehensive simulations of the two models demonstrate the potential application of this algorithm to the modeling of land-surface processes.

4 Conclusions and discussion

In this study, we use a two-directional freeze and thaw algorithm to simulate the FTFs which we incorporate into the thermal diffusion equation. A local adaptive variable-grid method is used to discretize the model. We perform sensitivity tests using different time steps and spatial resolutions, which demonstrate that the proposed method is stable and that the FTFs can be tracked continuously. We perform simulations of the FTFs and soil temperature at the Qinghai-Tibet Plateau D66 site from September 1, 1997 to September 22, 1998. The simulated values fit well with the data observations collected at the D66 site. The results show that the incorporated model performs much better in the soil temperature simulation than the original thermal diffusion equation, thereby demonstrating that this algorithm can potentially be used to the modeling of land-surface processes.

Acknowledgements The data used in this study are obtained from the Sino-Japanese international cooperation project “Global Energy and Water Balance Experiment —— Qinghai-Tibet Plateau Asian Monsoon Experiment” (GAME-Tibet).