Theory and Application of the Numerical Simulation in the Frozen Soil Problems

The freezing-thawing processes in soils are important components of terrestrial hydrology, which significantly influence energy and water exchanges between land surface and sub-surface. Long-term changes in frost and thaw depths are also an important indicator of climate change. A water-heat coupled movements model is established with frozen soil in this paper, which treats the freezing/thawing front as a moving interface governed by some Stefan problems with two free boundaries. The numerical simulation is conducted by using the modified finite difference method. The model is validated to compare its predictions with GEWEX Asian Monsoon Experiment(GAME)-Tibet observations at D66 site in Tibetan Plateau. The results show that the simulated soil temperature, soil water content and frost/thaw depth are in excellent agreement with the measured values. Finally, optimal error estimation for L∞ norm is derived on the model problem by using coordinate transformation method. The numerical simulation system is established on the basis of rigorous mathematics and mechanics, which successfully solved the important and difficult problems of environmental science.

to the climate system, which give delayed or rapid freezing/thaw.
In this paper, the seasonal depth changes in frozen and thawed soil is described to the coupled water and heat moving boundary problem on two free moving boundaries.The proposed model explicitly tracks the variation and dynamics of the active layer by computing the moving interfaces between the frozen zone and the unfrozen zone, and obtains the soil temperature, humidity and freeze/thaw depth simultaneously. According to the main physical processes in frozen soil, some hypotheses were proposed, including that: soil is an incompressible, homogeneous and isotropic medium; ice is immovable and only liquid water can move in frozen soil; water and heat movements in seasonally frozen soil only happen in the vertical direction; phase changes only take place on the phase-transition interfaces. Let (0, L) be a onedimensional vertical soil column from the fixed surface z = 0 at time t, and z = L is the bottom of the soil column, see Figure1. Phase-transition interfaces divide the soil profile into frozen and unfrozen regions. Thus, the whole soil column is supposed to be divided into three zones: a thawed layer from ground surface to the first phase-transition interface (namely thaw depth), a frozen layer from thaw depth to the second phase-transition interface (namely frost depth), and an unfrozen layer from frost depth to the soil bottom. Based on mass conservation, energy balance theory, and the temperature continuity in both phase change interfaces, the coupled soil water and heat problem in the freezing-thawing process can be described as follows: Where, T is soil temperature( • C); t, z are time(s) and depth(m)(positive downwards), respectively; ξ, ς are thaw depth and frost depth, respectively; c f , c u are thermal capacity(kJ · m −3 · C −1 ) for frozen and unfrozen soil, respectively. c f = (1−θ s )c s +θ i c i , c u = (1−θ s )c s +θc l , where c s , c i and c l are thermal capacity of dry soil, ice and liquid water, respectively, θ s is saturated water content; λ f , λ u are heat conductivity(kJ·m −3 ·C −1 ) for frozen and unfrozen soil, respectively(W·m −1 ·C −1 ).
For theory analysis of moving boundary problem (Stefan problem), see Crank (1984). Crank divides the numerical methods for solving moving boundary problems into 3 categories: boundary tracking method, boundary fixed method and region fixed method. To solve various moving boundary problems, some numerical methods have been developed in the last decades, such as finite difference method with coordinate transformation, finite element method, VOF method, integral equation method and the boundary element dual reciprocity method (see, e. g., Cannon(1984), Alexiades, V.V. & Solomon, A.D.(1993), Vermolen, F. & Vuik, K.(1998), Caldwell, J.& Savovic, S.(2002), Gupta (2003), Yuan (2011)). Some finite element methods (see Lunardini (1991)) have been presented for calculation and simulation of Stefan problems with various initial and boundary conditions (see, e. g., Asaithambin (1997),Černý, R.& Přikryl, P.(1999)) through continuous improvement. Segal et al. (1996) proposed an adaptive grid method with the use of the total time derivative. Wu (see Wu (2003)) presented wave front tracking method for the semi-infinite space. A large bibliography on the subject was given in (Tarzia, D.A.(2000)). However, most numerical methods and theoretical analysis focus on single phase or two phase Stefan problem in the semi-infinite space, and pay less attention to multi-phase problems or fixed region.
The paper is organized as follows. In Section 2, we adopt a finite difference scheme with fixed time step and local varied space step for the physical model. In Section 3, we do some experiments and compare the numerical results with GAME-Tibet observations at D66 site. The comparison shows that the simulated the soil temperatures, water content and soil frost/thaw depth have a similar trend to the measured values. Next, in Section 4, optimal rate of convergence of the model problem for the L ∞ norm is derived using coordinate transformation method. Finally, we conclude the paper with a summary in Section 5. The numerical simulation system is established on the basis of rigorous mathematics and mechanics, which is successful to solve the environmental science in this important and difficult problem.

Numerical Methods
Since soil water and heat coupled equations in the soil freezing and thawing process in the previous section are highly nonlinear and involve two unknown moving interfaces, it is very difficult to obtain analytical or semi -analytical solutions for the model. Consequently, the only feasible way to solve the problem is to compute approximate numerical solutions. We shall present finite difference scheme with local varied spatial grid nodes to obtain the numerical solutions. The frost depth and thaw depth of soil phase transition interfaces are taken as the grid nodes participating in the computation, and updated at each time step. Specifically, first obtain the soil temperature distribution from equations (1), (2), (4), (6). Then calculate the soil freezing and thawing depth by (8) and (9). Next, obtain soil moisture content from equations (3), (5), (6) and (10). Finally, update spatial grid nodes. The detail numerical scheme will be presented in the following.
Let τ = T N be time step size, {t n = nτ, n = 0, 1, · · · , N}, Ω τ = {t n |0 ≤ n ≤ N}. For [0, L], we adopt unequal interval mesh: 0 = z n 0 < z n 1 < · · · < z n J n = L, Ω n h = {z n j |0 ≤ j ≤ J n , 0 ≤ n ≤ N}. Moreover, let z = ξ and z = ς are just located some grid node j = s and j = r. Note that z = ξ and z = ς must be located in Ω n h at each time step. Because the possibility that calculated frost/thaw depth at time t = t n+1 coincides with the grid node in Ω n h is existence, so that J n is a positive integer varying with time. For brevity, we omit superscript n in the spatial grid nodes. Let h j = z j − z j−1 , j = 1, 2, · · · , J, h = max  Let {y n j = y(z j , t n )|0 ≤ j ≤ J, 0 ≤ n ≤ N} be the net function on the Ω h × Ω τ , define the quotients: Define L 2 norm and maximum norm of the function and its difference quotient: In this paper, we study the following difference scheme of equations (1)-(10): where ϑ and T are the approximate solutions of θ and T , respectively.
Initial and boundary conditions are discrete as follows: Therefore, we get the solvable three diagonal linear algebraic equations by the above difference scheme: where {x n+1 i |i = 1, 2, · · · , m} is the solutions of T or ϑ at time t n+1 on the grid nodes by (11)-(16); F n i = F i (ϑ n , T n , z i ), i = 1, 2, · · · , m. However, z = ξ and z = ς are always located in grid nodes at each time step, so it will be a set of variable dimension three diagonal linear algebraic equations with time. To implement the above finite difference methods on computers, we propose the following computation algorithm( the flow chart can be seen in figure3): Step 1. Input initial data: choose soil temperature, soil water content and boundary conditions g 1 (z), g 2 (z), f 1 (t), f 2 (t) on the initial time. Set n = 0, let ξ 0 = ς 0 = 0, then s = r = 0 at time t = t 0 ; Step 2. Calculation soil temperature: obtain T n+1 j from scheme (12), (13) according to (17); Step 3. Computation frost/thaw depth: obtain ς n+1 and ξ n+1 from scheme (14), (15); Step 4. Calculation soil water content: obtain ϑ n+1 j from scheme (11) according to (17); Step 5. Update space grid nodes: update Ω h = {z j | j = 0, 1, · · · , J} based on ξ n+1 , ς n+1 . First, delete ξ n , ς n at t = t n from Ω h . Next, add ξ n+1 and ς n+1 calculated by step3 into Ω h . Then space grid nodes at time t = t n+1 can be taken; Step 6. Let n = n + 1 and return to step2. Continue the time recycle until n = N.

Ideal Experiment
In this section, set L = 2.5m, which is equidistant divided into 20 layers on the beginning. Assume soil moisture is remained unchanged at 0.06 in the simulation process. The upper boundary condition for temperature is taken as the following sinusoidal periodic boundary conditions: where T 0 = 0 • C is mean annual ground surface temperature(GST); G t = 0.02 • C/a is the rising rate of GST; A 0 = 13 • C is annual amplitude of GST, i.e., half of annual GST differences; ω = 2π 8760 . The zero heat flux boundary condition at the bottom of the soil column are used. Let , and S=0. Figure4 shows the simulated results. Figure 4 shows soil surface temperature and simulated frost/thaw depth varying with time. As you can see from Figure  4, there is obviously sinusoidal/cosine variation of soil frost and thaw depths with the sinusoidal variation of land surface temperature on one year timescale. At the same time, frost/thaw depth also has obvious inter-annual variability with the annual period change of land surface temperature on the interannual timescale. The phase difference between the period of the land surface temperature and that of frost depth is 8 days, and the phase difference between that of the surface temperature and that of the thaw depth is 4 days. Meanwhile, the minimum frost depth and the maximum thaw depth become deeper with surface temperature increasing. It shows that the long-term changes of the frost and thaw depths are affected by the interannual change of surface temperature. Therefore, the model is reliable.

Numerical Simulations With Observation Data
In this section, we use the data from the Japan international cooperation project "global energy and water balance experiment -the Asian monsoon in Qinghai Tibet Plateau" (GAME-Tibet) to validate the numerical model. D66 site is located in the Qinghai-Tibet highway with an elevation of 4610 m, and the landscape at this site is short-grazed grassland. We take the upper boundary conditions as the observations of surface soil temperature and humidity in the simulation process. And the zero heat flux and residual water content are used as the lower boundary condition. Meanwhile, the observations on the beginning are taken as the initial conditions. The model was tested using the data from 1997.8.19 to 1998.8.31 at the D66 site. The values of the parameters are defined as follows( see Li and Koike (2003)).
Figure5 shows the simulated soil frost depth and thaw depth from 1997.8 to 1998.8 at D66. The simulation differences are mainly focused on the freezing period in winter. It can be seen that D66 site begins to freeze in October and melt in April on the surface. The period of the freezing-thawing process is about 6 months, and completes this process in mid June. As can be seen, the freezing process takes about two months and freezing depth can reach about 200cm in late November according to the simulated frost depth. After this, the freezing velocity is obviously slowed down, and the growth rate of freezing depth is no more than 50cm. Soil frost depth can reach about 250cm in early February after 4 months. In April, shallow soil has began to melt and continuously passed down despite the deep soil continues to freeze. Model simulation results can reflect the process well.
Figure6 shows the simulated and observed soil temperature at 4cm, 20cm and 100cm respectively. We can see that soil temperatures complete a cycle in a year period. In 1998 January, the simulated soil temperature has dropped to the lowest point (252.3K, 254.55K, 263.28K), and then began to warm. Beginning in February, shallow soil temperature has started to rise, and reached the highest value in a year (290K, 288K, and 281K) in late June. From the beginning of July soil temperature has decreased. To the simulated soil temperature from August to December in the cooling process, the cooling rate of each month is large. The monthly average temperature range is more than 6K. From February to June, the monthly average temperatures are greatly warmer in the range of about 5K. The monthly variations of soil temperature of each month from June to August become smaller. These results are consistent with the observations. At the same time, soil temperatures show sinusoidal change in figure 6(e) obviously, which is consistent with the annual variation of solar radiation. The long-term changes in the upper soil temperature is mainly influenced by the interannual variation of solar radiation.
Figure7 shows the differences between the simulation and observation of soil temperature at site D66. It can be seen that the simulation results and observations are still relatively close. In the shallow layer of soil (4cm), most differences are not more than 3 • C besides few simulation results. Moreover, temperature differences are between −2 • C and 2 • C in 20cm and 100cm. The average errors of simulation results in each layer are -0.0015, 1.2186e-004 and -0.0015, respectively. Since 1998 May the surface in 4cm has started to thaw and water content began to increase. As you can see from figure8, simulated error is mostly in 6-7 month, and most violently in the 4cm depth. The mean errors of 4cm, 20cm and 100cm are 0.0538, 0.0131 and 0.0027, respectively. At the same time, it can be seen that there is a negative correlation between the soil temperature at 4cm and water content changed with time from figure8(d).

Numerical Analysis
Due to the complexity of the practical problems, we only consider the numerical analysis of the model problem for equations (1) -(10) in this section. Thus, the model problem is equations (1)- (7), and the corresponding difference scheme is (11) -(13), (16). We assume that the conditions (I)-(V) are as follows: (I) ξ(t) and ζ(t) are the known functions depended on time t, ξ(t), where L * , ξ * , ξ * * , ς * are all positive constants; (II) λ u , λ f , c u , c f , D and K only depend on z, and do not depend on the unknown function θ, T and time t; (III) λ u and λ f have positive upper and lower bounds λ * and λ * , that is, 0 < λ * ≤ λ u , λ f ≤ λ * ; (IV) c u and c f have positive upper and lower bounds c * and c * , that is, 0 < c * ≤ c u , c f ≤ c * ; (V) D has positive lower bounds d * , that is, D ≥ d * > 0;

Coordinate Transformation
In order to analyze the convergence of the above scheme, let us introduce the change of variable: (1)- (3), we can write the following form: Accordingly, the initial and boundary conditions of eqs.(20), (21) are the followings: Let {x i , i = 0, 1, 2, · · · , J} are space nodes, h i = x i − x i−1 is space step size. x = 1 and x = 2 are just located on the node x s and x r . Let u n j = u(x j , t n ) and φ n j = φ(x j , t n ) are the approximate solutions to U n j = U(x j , t n ) and θ n j = θ(x j , t n ), respectively; ς n = ς(t n ), ξ n = ξ(t n ).
Theorem1: Suppose the assumption conditions (I)-(V) hold. Let U and θ be the solutions of the model problem (20)-(21) with the initial and boundary conditions (22), u and φ solve the finite different scheme (23)-(25). Then for U, θ ∈ W 2,∞ (L ∞ ) ∩ L ∞ (H 2 ), when τ and h tend to zero at the same time, we have where M * is a constant that is not dependent on h and τ, but dependent on the functions U and θ and its derivatives.
Furthermore, we get the following second order L 2 norm estimation by the detailed analysis and deduction if uniform mesh generation used.
Theorem2: Under the conditions of Theorem1, the errors in algorithm (23)-(25) can be estimated by: where f 2 x, j h j , M * * is a constant that is not dependent on h and τ, but dependent on the functions U and θ and its derivatives.

Conclusions
In this paper, water and heat coupled process in the frozen soil is attributed to Stefan problems with two moving boundary conditions. The finite difference scheme is presented with fixed time step and variable space step. It is shown that the model is effective and reasonable, the calculation process is simple and physical meaning is clear by comparison numerical simulation results with observations. Finally, the algorithm error is analyzed with coordinate transformation, and the optimal L ∞ error estimate is obtained. This thesis will be further studied in high dimensional problems on the basis of the above work.