A Model for Coupled Moisture and Heat Flow in Unsaturated Soil S. Kim Manager, Ph.D., Dept. of Project Development, Korea Land Corporation, 217 Jeongja-Dong Pundang-Gu Sungnam Kyunggi 463-755, Korea and A. G. Heydinger Professor, Department of Civil Engineering, University of Toledo, Toledo, OH, USA

Abstract

Soils within and below pavement systems are designed to be in an unsaturated state. Therefore, moisture and temperature distributions in unsaturated soils are affected by climatic factors. Hydraulic gradients due to rainfall cause moisture flow and sensible heat transport, and thermal gradients due to temperature changes induce not only heat flow but also moisture flow. Freezing temperatures significantly reduce the permeability of soil but increase moisture flow caused by hydraulic gradients due to ice lens formation in the frozen soil that act as sinks.

Heat flow and moisture flow have been recognized as coupled processes with complex interactions between them. This paper presents an analysis model by the finite element method for the coupled heat flow and moisture flow in unsaturated soils. The model is applicable to both one-dimensional and two-dimensional problems. The model can be used to predict not only the change of temperature and water content, but also frost heave with time. The model is tested through comparisons with results by other computational models and by experimental results.

Keywords: coupled flow analysis, finite element method; permeability function; soil-water characteristic function; unsaturated soil.

Introduction

The distribution of moisture and temperature within pavement systems vary seasonally depending on environmental conditions. Since moisture and temperature are important factors for pavement design, it is important to evaluate the seasonal variations of moisture and temperature in pavement systems. Moisture and heat flow are dependent processes. Therefore, mathematical models that fully couple moisture and heat flow are required to model moisture and temperature changes within unsaturated soils. Models that predict moisture or temperature changes independently do not adequately represent physical phenomena that occur. However, they are useful for verification purposes as was done for this model. One-dimensional mathematical models that include coupled moisture and heat flow in the vertical direction only are not adequate for modeling moisture and temperature variations under pavements if there is significant horizontal variation of moisture and temperature. The present model described in this paper allows for a fully coupled analysis of one-dimensional or two-dimensional flow in unsaturated soils with frozen or unfrozen conditions.

Two empirical equations that give the material dependent relationships relating soil moisture and permeability to soil suction are required to solve the governing equations of the coupled moisture flow. One is the relationship between water content and suction for the soil, the soil-water characteristic equation. The other is the relationship between permeability and suction, the permeability function. Researchers have proposed many equations for these relationships. The present model uses two sets of equations that are the soil-water characteristic and permeability functions. One set consists of the equations by Gardner (1958) and Guymon et al. (1993) and the other set consists of the equations by Fredlund and Xing (1994) and Fredlund et al. (1994).

Literature review

Although the movement of moisture and heat in unsaturated soils has been of interest to researchers since the early 1900’s, the theoretical development for the movement of moisture and heat has been available for a few decades. Attempts to analyze the coupled heat and moisture flow using mathematical methods have been presented by Harlan (1973), Guymon and Luthin (1974), Jame and Norum (1980), Hromadka et al. (1981 a and b), Guymon et al. (1993) and Newman and Wilson (1997). The models are similar; they are based on some form of the equations given by Harlan (1973).

There are differences in the applications of the numerical equations and soil properties. For example, Harlan (1973), Jame and Norum (1980), and Newman and Wilson (1997) neglected the convective (or advective) term of the heat equation but Guymon and Luthin (1974), and Guymon et al. (1981, 1993) considered the term. Harlan (1973) used the finite difference method to model the one-dimensional coupled heat and moisture flow. He assumed that the soil-water characteristic curve and the suction versus hydraulic conductivity relationship would be same in the unfrozen and the frozen soils. Jame and Norum (1980) used the finite difference methods like Harlan (1973) but they noted that there was some difference between the numerical results based on Harlan (1973) and the results obtained from their measurements. They thought the difference was caused by the ice accumulation that disrupted the flow paths and hence reduced the flow rate. Jame and Norum (1980) introduced an impedance factor to account for the reduced flow in the frozen soil.

A one-dimensional coupled heat and moisture flow modeling using the finite element method was developed by Guymon and Luthin (1974), Guymon et al. (1981, 1981), Newman and Wilson (1997), and Larson and Dempsey (1997). Guymon et al. (1993) used Gardner’s (1958) equation to establish the soil-water characteristic curve and the suction versus hydraulic conductivity relationship. They also thought ice might be partly blocking soil pores and reduce hydraulic conductivity, and thus used an E parameter similar to the impedance factor proposed by Jame and Norum (1980). Newman and Wilson (1997) used the soil-water characteristic curve equation presented by Fredlund and Xing (1994) and the suction versus hydraulic conductivity relationship presented by Fredlund et al. (1994). They did not use an impedance factor to consider the decrease of hydraulic conductivity due to ice. Newman and Wilson (1997) thought that the permeability function by Fredlund and Xing (1994) is an improvement to permeability functions used by earlier researchers. They showed that there was good agreement between their model predictions and experimental results.

Larson and Dempsey (1997) published a model called the Enhanced Integrated Climatic Model (EICM) Version 2.0 that is an update to the Integrated Climatic Model (ICM) developed by Lytton et al. (1993). The ICM is a one-dimensional coupled heat and moisture flow model to analyze pavement soil systems subject to climatic effects. The ICM uses the relationships that were used by Guymon et al. (1993) for the soil-water characteristic curve and permeability functions. In addition, Lytton et al. (1993) used a method presented by Washburn (1924) for estimating the suction in the frozen zone to solve the problem for the decrease of hydraulic conductivity instead of using an E factor.

The EICM allows for the use of the equations by Gardner (1958) and Guymond et al. (1993).

The solutions described above are one-dimensional numerical models for the coupled heat and moisture flow. A numerical model described by Hromadka et al. (1981b) presented an analysis procedure for two-dimensional moisture flow in an unsaturated soil using the finite element method where the solution domain was divided into triangular elements. Two and three-dimensional models for uncoupled moisture (SEEP/W) and temperature (TEMP/W) flow are available as commercial software packages (Geo-Slope International Ltd).

Moisture and heat flow equations

Two systems of equations are required to describe mathematically the coupled heat and moisture flow. One is the moisture flow equation that is an application of the principle of mass conservation and the other is the heat flow equation that is an application of the principle of energy conservation.

The two-dimensional moisture flow equation in an unsaturated soil is

 (1)

where kx(y) and ky(y)are the hydraulic conductivities, y is the suction, c(y) is specific moisture capacity, ri is the density of ice, rw is the density of liquid water, and qi is the volumetric ice content. The specific moisture capacity is defined as

 (2)

To solve equation (1) the specific moisture capacity cw and the hydraulic conductivities kx(y) and ky(y) must be known. The specific moisture capacity can be determined from the soil-water characteristic curve and the hydraulic conductivities can be determined from the permeability function. In the model, two sets of equations were used. One is soil-water characteristic curve by Gardner (1958) and permeability function by Guymon et al. (1993). These equations are relatively simple and thus have been used by many researchers. The other set is soil-water characteristic curve by Fredlund and Xing (1994) and permeability function by Fredlund et al. (1994). These equations are shown in the appendix.

The two-dimensional heat flow equation is

 (3)

where is KT thermal conductivity, T is temperature, vx and vy are the velocities of the water flow in x and y directions, Cw is the heat capacity of water, Cm is the volumetric heat capacity, L is the latent heat of fusion of water which is the thermal energy required for a change of phase, ri and rw are the densities of the ice and water, and qi is the volumetric ice content. The thermal parameters in equation (3) can be computed from DeVries (1966) relationship:

 (4)

 (5)

where qu is unsaturated volumetric water content, q0 is porosity of soil, Cw, Ci and Cs are volumetric heat capacities of water, ice and soil, and Kw, Ki and Ks are thermal conductivities of water, ice and soil. The velocity of the water flow is calculated from Darcy’s law:

 (6)

The last terms of equations (1) and (3) are sink or source terms due to the change of ice phase. These terms are approximated by an isothermal phase change process presented by Hromadka et al. (1981a). Then, dropping off the last terms of equations (1) and (3) become

 (7)

 (8)

Finite Element Method Solution

The following equations are the numerical equations for equations (7) and (8) which were derived by the finite element method.

 (9)

 (10)

In equation (9), [C (y)] is called the capacitance matrix and [K(y)] is called the conductance matrix. w is a value depending on the numerical method: It is 0 for the forward difference method, ½ for the central difference, and 1 for the backward difference method. In this study the backward difference method was used. Also, {F} is the specified rate of the flow. In equation (10), [C] is the capacitance matrix and [D] is the conduction-advection matrix.

Programming FEM Solution

The matrix equations derived for the finite element method are shown in equations (9) and (10) (Kim, 2002). To program the FEM solution a programming language that deals efficiently with matrix equations is needed. MATLAB was used for this study (Hanselman and Littlefield, 1997). The analysis procedure in the program is: first, initial and boundary conditions and soil parameters are set up; second, moisture flow is analyzed; and third, heat flow is analyzed. The results analyzed by moisture and heat flow are then corrected by an isothermal phase change approximation (Hromadka et al., 1981a). Finally the corrected results are saved. This procedure continues until the end of each time step.

Initial and Boundary Conditions

For the program, initial and boundary conditions are specified as input data. The initial conditions specified are water content, temperature and ice content. For the boundary conditions, moisture and temperature values are specified. The moisture values specified are suction, volumetric water content or water flux. Total head boundary conditions were assumed for the top and bottom boundaries. In case of two-dimensional problem, both sides of pavement section were assumed as no flow boundaries.

Structure of model

Figure 1 shows the structure of the model. The model consists of four major components. They are command window, input, output, and the component for analysis of coupled heat and moisture flow.

Figure 1. Structure of Model

Model verification

The model presented here was verified through comparisons with analysis results by other models and with measured data. For the comparisons with the analysis results, it was necessary to verify the model with models that do not allow coupled moisture and heat flow. The assumption is that, if the present model is verified through comparisons with uncoupled models, the model may be used for coupled analysis with confidence. Three model verification examples are shown in this paper.

Example 1

Philip (1957) solved infiltration problems in Yolo light clay using a quasi-analytical method. The soil-water characteristic curve and permeability function used in Philip’s solution were expressed by Haverkamp et al. (1977) as the following equations:

 (11)

 (12)

where Ks = 0.04428 cm/hr

A problem was solved with the following boundary and initial conditions:

 t = 0   0 >=z >= -50cm   q =0.2376 (or y = -600) t > 0   z=0   q = 0.495 (or y = 0) t > 0   z = -50   q = 0.2376 (or y = -600)

This problem is a one-dimensional problem. The problem domain was divided into 40 elements where the top 30 elements have 1-cm length and the bottom 10 elements have 2-cm length. The time step size is 0.2 hour. Figure 2 shows the results by the present model and Philip’s solution.

Figure 2. Comparison Between the Results by the Present
Numerical Model and Philip’s Analytical Solution

Example 2

Seasonal variations of moisture and temperature within pavements currently being monitored with the Seasonal Monitoring Program (SMP) instrumentation prescribed by Strategic Highway Research Program (SHRP) protocols which is included as a part of the Long Term Pavement Performance (LTPP) testing (FHWA, 1994; Heydinger and Randolph, 1998). Figure 3 shows a pavement section of U.S. Route 23 in Delaware County, Ohio where instrumentation sensors were included. The pavement section consists of 27.9 cm Portland cement concrete (PCC), 15.2 cm graded aggregate base (DGAB) and subgrade soil.

Figure 3. A Pavement Section of U.S. Route 23 In Delaware County, Ohio
(From Ohio University, 1994)

The results by monitoring and by the model presented here are shown in Figure 4. The results show relatively good agreement as shown in the figures. At a depth between 1.2 and 1.4 meters shown in Figure 4 (a) the measured water contents are less than the simulated water contents. It is considered that there is some inhomogeneous subgrade soil, though it was assumed that the subgrade soil is homogeneous for the numerical analysis.

(a) Water Content Change Analysis In The Subgrade Soil

(b) Temperature Change Analysis

Figure 4. Comparison Between The Results
By Measurement and New Model

Example 3

The model presented here has been developed for two-dimensional analysis as well as one-dimensional analysis for coupled heat and moisture flow. The verification of the two-dimensional model was performed using by GEO-SLOPE which is a software developed by GEO-SLOPE International Ltd. (GEO-SLOPE, 1991-1998). GEO-SLOPE OFFICE consists of computer programs called CTRAN/W, SEEP/W, SIGMA/W,

SLOPE/W, and TEMP/W. SEEP/W and TEMP/W can be used for the verification. SEEP/W is a seepage analysis program that models both saturated and unsaturated flow. TEMP/W is a program that models the thermal changes in the ground.

Figure 5 shows the problem domain for numerical analyses that consists of two layer soils and a drainage pipe. The problem domain was divided into 86 elements and has 108 nodes. Figure 6 shows the analysis results for water content changes with time by the two models. The figures on the left side are the analysis results by SEEP/W and the figures on the right side are the analysis results by the model presented here. Figure 7 shows the analysis results for temperature changes with time by the two models. The figures on the left side are the analysis results by TEMP/W and the figures on the right side are the analysis results by the model presented here.

Figure 5. Problem Domain

Figure 6. Comparison of Water Content Changes by Two Numerical Models

Figure 7. Comparison of Temperature Changes by Two Numerical Models

Summary and Conclusions

An analysis model for the coupled heat and moisture flow was presented. The model presented here can be used for both one-dimensional and two-dimensional problems. It is considered that two-dimensional analysis may be needed more than one-dimensional analysis in practice, because water content and temperature might change horizontally or vertically.

Two equations for the relationships between soil parameters are used for moisture flow. They are the soil-water characteristic relationship and the permeability functions. Two sets of equations were used by the present model One is the soil-water characteristic relationship by Gardner (1958) and the permeability function by Guymon et al. (1993). The other set is soil-water characteristic relationship by Fredlund and Xing (1994) and permeability function by Fredlund et al. (1994).

The model presented here was developed to analyze seasonal variations of water content and temperature as well as frost heave prediction. In the present work, the verification for the analysis of water content and temperature variation within pavement has been tested. The results by the model show good agreement with the results by others. There is a need to verify the model for frost heave in instrumented pavements.

References

1. DeVries, D. A. (1966) “Thermal Properties of Soils.” In Physics of Plant Environment (W. E. Van Wijk, Ed.), Amsterdam: North-Holland Publishing Co., 210-235.
2. Federal Highway Administration (1994), “LTPP Seasonal Monitoring Program: Instrumentation Installation and Data Collection Guidelines.” FHWA-RD-94-110, Research and Development, Turner-Fairbank Highway Research Center, 6300 Georgetown Pike, McLean, VA 22101-2296.
3. Fredlund, D.G. and A. Xing (1994) “Equation For The Soil-Water Characteristic Curve.” Canadian Geotechnical Journal, 31, 521-532.
4. Fredlund, D.G., A. Xing and S. Huang (1994) “Predicting The Permeability Function For Unsaturated Soils Using The Soil-Water Characteristic Curve.” Canadian Geotechnical Journal, 31, 533-546.
5. Gardner, W.R. (1958) “Some Steady State Solutions of the Unsaturated Moisture Flow Equation with Application to Evaporation From a Water-Table.” Soil Science, 85, 228-232.
7. Guymon, G. L., R. L. Berg and T. V. Hromadka (1993) “Mathematical Model of Frost Heave And Thaw Settlement In Pavements.” CRREL Report 93-2.
8. Guymon, G. L., R. L Berg, T. C. Johnson and T.V. Hromadka (1981) “ Results From A Mathematical Model of Frost Heave.” Transportation Research Report 809, 2-6.
9. Guymon, G. L. and J. M. Luthin (1974) “ A Coupled Heat And Moisture Transport Model For Arctic Soils.” Water Resources Research, Vol. 10, No. 5, 995-1001.
10. Hanselman, D. C. and B. Littlefield (1997), “The Student Edition of MATLAB:Version 5, User’s Guide, The MathWorks, Inc, Upper Saddle River, NJ, 285-294.
11. Harlan, R. L. (1973) “ Analysis of Coupled Heat-Fluid Transport In Partially Frozen Soil.” Water Resources Research, Vol. 9, No. 5, 1314-1323.
12. Haverkamp, R., M. Vauclin, J. Touma, P.J. Wierenga, and G. Vachaud (1977) “ A Comparison of Numerical Simulation Models For One-Dimensional Infiltration.” Soil Science Society of America Journal, Vol. 41, No. 2, 285-294.
13. Hromadka, T. V., G. L. Guymon and R.L Berg (1981a) “Some Approaches To Modeling Phase Change In Freezing Soils.” Cold Regions Science and Technicalogy, 4:137-145.
14. Hromadka, T. V., G. L. Guymon and G.C. Pardoen (1981b) “ Nodal Domain Integration Model of Unsaturated Two-Dimensional Soil-Water Flow: Development.” Water Resources Research, Vol. 17, No. 5, 1425-1430.
15. Heydinger, A. G. and B. W. Randolph (1998) “Seasonal Instrumentation of SHRP Pavement – University of Toledo.” Report No. FHWA/OH-98/012.
16. Jame, Y. W. and D. I. Norum (1980) “Heat And Mass Transfer In A Freezing Unsaturated Porous Medium.” Water Resources Research, Vol. 16, No. 4, 811-819.
17. Kim, S. (2002), “Coupled Heat and Moisture Flow Analysis in Unsaturated Soil.” Ph. D. Dissertation, Department of Civil Engineering, University of Toledo, Toledo, OH.
18. Larson, G. and B. J. Dempsey (1997), “Enhanced Integrated Climatic Model Version 2.0.” DTFA MN/DOT 72114.
19. Lytton, R. L., C. H. Pufahl, H. S. Michalak, H. S. Liang and B. J. Dempsey (1993), “An Integrated Model of the Climatic Effects on Pavements.” FHWA-RD-90-033.
20. Newman, G. P. and G. W. Wilson (1997) “Heat And Mass Transfer In Unsaturated Soils During Freezing.” Canadian Geotechnical Journal, 34: 63-70.
21. Ohio University (1994), “Instrumentation Plan for the Ohio SPS Test Pavement (DEL-23-17.48).” Center for Geotechnical and Environmental Research, Civil Engineering Department.
22. Philip, J. R. (1957) “The Theory of Infiltration Equation And Its Solution, Soil Science, 83(5), 345-357.
23. Washburn, E. W. (1924), “The Vapor Pressure of Ice and Water Below the Freezing Point.” Monthly Weather Review, National Research Council, Vol. 52, 488-490.

Appendix

Soil-water characteristic curve by Gardner (1958):

where qu is the unsaturated water content, qs is the volumetric water content at saturation, Aw and a are best fit parameters, and y is the suction.

Soil-water characteristic curve by Fredlund and Xing (1994):

where C(y) is a correcting factor, e is exponent value, and a, n and m are best fit parameters.

Suction versus permeability relationship by Guymon et al. (1993):

where hh is the unsaturated hydraulic conductivity, ks is the saturated hydraulic conductivity, and Ak and b are best fit parameters.

Suction versus permeability relationship by Fredlund et al. (1994): Fredlund et al. presented an equation for the relative coefficient of permeability to predict the permeability of unsaturated soil using the soil-water characteristic curve equation by Fredlund and Xing (1994).

where kr is the relative coefficient of permeability, b = ln(1000000), is the derivative of equation by Fredlund and Xing (1994) with respect to y, and yaev is the air entry value. kr is the ratio of the unsaturated permeability to the saturated permeability in a soil.