Transport coefficients and pressure conditions for growth of ice lens in frozen soil

In this paper, the transport of sub-cooled water across a partially frozen soil matrix (frozen fringe) caused by a temperature difference over the fringe, is described using non-equilibrium thermodynamics. A set of coupled transport equations of heat and mass is presented; implying that, in the frozen fringe, both driving forces of pressure and temperature gradients simultaneously contribute to transport of water and heat. The temperature-gradient-induced water flow is the main source of frost heave phenomenon that feeds the growing ice lens. It is shown that three measurable transport coefficients are adequate to model the process; permeability (also called hydraulic conductivity), thermal conductivity and a cross coupling coefficient that may be named thermodynamic frost heave coefficient. Thus, no ad hoc parameterizations are required. The definition and experimental determination of the transport coefficients are extensively discussed in the paper. The maximum pressure that is needed to stop the growth of an ice lens, called the maximum frost heave pressure, is predicted by the proposed model. Numerical results for corresponding temperature and pressure profiles are computed using available data sets from the literature. Frost heave rates are also computed and compared with the experimental results, and reasonable agreement is achieved.


Introduction
Frost heave, i.e. the transport of sub-cooled water to a growing ice lens, occurs when three conditions coincide: the temperature is below the normal freezing point of water, the sub-cooled water is connected to a water reservoir, and the soil is susceptible to formation of ice lens. The growth, which often takes place under confined conditions, can lead to a substantial build-up of pressure, or ground self-organization [11]. A maximum value of 1.14 MPa/K was reported by Takashi et al. [26] and Førland and Kjelstrup Ratkje [8]. This pressure can result in service life reduction, and sometimes failure of the overlaying infrastructures. Criteria and properties that characterize the phenomenon are therefore of interest.
In a series of papers, Konrad and Morgenstern [17,19,18] elaborated on a physical model for frost heave. They divided the soil into three layers; (1) a layer of unfrozen soil, (2) a partially frozen layer, called the frozen fringe, and (3) a practically completely frozen one. These layers are illustrated in Fig. 1. The temperature and pressure at the top of the unfrozen layer are T 0 (0°C or 273.15 K) and p 0 , respectively. Inside the frozen fringe (layer 2) T s \T\T 0 . At the base of the ice lens the temperature is T s , the ice segregation temperature. Layer (3) above the ice has T\T s and is completely frozen. Similar terminology is adopted in this paper (Fig. 1).
Assuming continuity of the subcooled liquid water phase, the liquid pressure varies from p 0 at the bottom of the frozen fringe to p l at the bottom of the ice lens, as shown in Fig. 1. The frozen fringe thickness is determined by the freezing temperature of the subcooled confined water. In coarser soils, the position of the ice lens can be close to the freezing front (a few mm), while in finer soils, it can be many millimeters away from the freezing front.
Throughout this work, the common Darcy's and Fourier's laws are adopted for layers 1 and 3. While, in the frozen fringe, layer 2, the mass flux carries heat. This implies coupled transport of heat and mass, and thus the simple Darcy's and Fourier's laws are not valid for this layer. By coupling, we mean specifically that the mass flux might also be driven by a thermal force, and likewise, that the heat flux will obtain a contribution from the pressure gradient. This coupled transport phenomenon was already proven by experimental evidence [17,24,16].
The proportionality coefficient between the water flux and the temperature gradient was measured and introduced as the segregation potential (SP), by Konrad and Morgenstern [18]. However, they used this formula in the presence of a pressure gradient. While, as shown, for a well-defined measurement, it is critical to specify a constant pressure along the frozen fringe. To implement this modification, we propose that the thermodynamic frost heave coefficient, s, replaces SP.
Non-equilibrium thermodynamics theory was already applied to frost heave problem in 1967. According to the review of Black and Hardenberg [3], Hoekstra and Miller [12] began to incorporate concepts of irreversible thermodynamics. The review also cited Groenevelt and Bolt [9], Kay and Groenevelt [13], and Groenevelt and Kay [10] for use of a similar approach. The rigid-ice model of Miller [21] described the transport of heat and mass through a frozen soil using similar equations from irreversible thermodynamics. The theory was also applied by Derjaguin and Churaev [6], and Førland and Kjelstrup Ratkje [8] to describe the phenomenon. More recently, a similar formulation, but from a different approach, was derived by Rempel [25], and Wettlaufer and Worster [29]. All authors predicted that a maximum pressure exists, and such a pressure has been measured by Takashi et al. [26]. A possibility to improve on the description of frost heave is nevertheless possible using the discrete theory of nonequilibrium thermodynamics for heterogeneous systems [14]. This is important because of the discrete nature of the process, and the need to find an equation set suitable for modelling on larger scales. The system pictured in Fig. 1 is discrete, and we need to deal with the overall conditions, the finite thickness of the frozen fringe, as well as the local conditions at the point of the ice lens.
The purpose of this article is thus to recapitulate and extend the earlier description of frost heave in terms of non-equilibrium thermodynamics. The aim is to achieve a formulation that can be used to obtain a necessary and sufficient set of transport coefficients from Onsager reciprocal relations, a set that later can be used in frost heave simulations. This work will therefore focus on specifying Onsager coefficients that characterize frost heave phenomena, and their laboratory measurements. Once the properties are known, mostly from measurements at steady state, they can be used to simulate the dynamic process. Such analyses will follow in work to come. Numerical examples will be given as far as possible, using properties of the well-studied Devon silt [18,1], and Xuzhou silty clay [33]. To have well-defined experiments, we continue to use the same three-layer characterization as Konrad and Morgenstern [17], and limit ourselves to water-saturated soils. The outline of the paper is as follows: The general coupled transport equations for the finite frozen fringe are defined in Sect. 2. A detailed description of the experimental determination of the coefficients follows in Sect. 3. A practical set of the transport equation, tailored for frost heave phenomena, is introduced in Sect. 4, and is followed by discussion and numerical examples illustrated with data from literature.

Coupled transport equations
The original layer construction and terminology of Konrad and Morgenstern [17] are adopted in this paper. Measurements are frequently done at steady state. Temperatures are in practice measured with an accuracy of 0.1 K and will be given accordingly. As explained in connection with Fig. 1, transport of sub-cooled water takes place from the unfrozen layer (at 273.2 K and some liquid pressure p 0 ) to an ice lens at a finite distance away. The ice lens is allowed to grow above the frozen fringe, at temperature T s \T 0 and Fig. 1 Schematic illustration of frost heave in layer 3 after transport of heat and water through layer 2. Layer 1 is connected to a water reservoir liquid pressure p l . The temperature across the frozen fringe may cover the range from 273.2 to 269.2 K for fine clay, or from 273.2 to 273.0 K, for coarser materials. For the Devon silt, the range 273.2-272.5 K was estimated by Konrad and Morgenstern [17].
Considering a volume element along the x-axis in Fig. 1, transport of heat and water takes place from the unfrozen layer towards the ice lens. It is assumed that the sample is fully consolidated prior to the steady state condition, and thus, there is no change in effective stress during the analysis. In layers 1 and 3, the simple laws of Darcy's and Fourier's are adopted. The coupling between transport phenomena takes place in layer 2. It is described by a symmetric flux-force matrix, standard of nonequilibrium thermodynamics [14]. We refer to the literature [14] for the basic principles and the derivation of the entropy production.
Following Kjelstrup and Bedeaux [14], the flux of heat _ q (in J m À2 s À1 ) and mass _ w (in kg m À2 s À1 ) out of layer 2 into layer 3 are: The fluxes are considered to be one-dimensional along the x-coordinate, and the difference D refers to the frozen fringe, which is a discrete system. The fringe includes water-ice phase transitions at both boundaries. We therefore always integrate from T 0 = 273.2 K and p 0 . The end of the fringe is the position where the ice lens is formed. The thermal driving force, DðT À1 Þ Dx (in K À1 m À1 ), is the main driving force for the heat flux _ q. Theory prescribes that this flux is taken at layer 3, when the driving force, The temperature in these equations is given in K. The chemical driving force has only one contribution in the absence of solutes; namely from the pressure change in the liquid water column across layer 2: Here V w (in m 3 kg À1 ) is the specific volume of water, while p is the fluid pressure. The difference refers to the pressure difference below the ice lens, and the pressure at the table of water. According to Kjelstrup and Bedeaux [14], _ q refers to the heat flux out of the layer 2 into layer 3, when T 0 is used in the definition of the chemical driving force. Derjaguin and Churaev [6] who used the same form, explained that the disjoining pressure did not have an impact on the form.
The coefficients L ij are generalized conductivity coefficients, to be related to more common coefficients in Sect. 3. The subscripts q and w refer to the phenomena they describe, the heat and water fluxes, respectively. The coupling coefficients L qw and L wq are special for irreversible thermodynamics. The coefficient matrix is symmetric according to Onsager, L qw ¼ L wq , but only if the fluxes and forces are written in the special form dictated by the entropy production [14]. Only three coefficients are independent, thanks to the matrix symmetry. In Sect. 3, the independent coefficients will be defined and be related to experiments. The coefficients obey D ¼ L qq L ww À L qw L wq ! 0 according to the second law of thermodynamics [14], and this relation can be used to check if coefficients are consistent with the second law.
The transport coefficients in Eqs. (1) and (2) depend on the state of the soil. Saturated frozen soils are in general characterized by their pore size distribution, mineralogy of the grains, unfrozen water saturation,Ŝ w , (alternatively, ice saturation,Ŝ i ) and temperature T. The presence of possible solutes in the water may also play a role. For a given soil composition, the unfrozen water saturation can be expressed by a freezing characteristic curve. Konrad and Duquennoi [16] gave a relation for Devon silt. In their model, the unfrozen water saturation was written as a logarithmic function of temperature Azmatch [1] reviewed the procedures for its determination and concluded that the unfrozen water saturation as well as permeability of frozen soils can be taken as a unique function of temperature (ignoring the hysteresis effect). When S w ¼ f ðTÞ, the transport coefficients will also be functions ofŜ w (or T) alone; i.e. L ij ¼ L ij ðTÞ.

Experimental determination of transport coefficients
Equations (1) and (2) will be used to define experiments. We proceed to explain a set of experiments that determines the Onsager coefficients, before we conclude with a more practical set of equations, containing more familiar transport coefficients. It should be noted that the samples to be measured in the following experiments are fully saturated with water before the experiment starts, and consolidated. The apparatuses are schematically shown in Figs. 2, 3 and 4.

Hydraulic conductivity (permeability)
As mentioned earlier, ice lenses usually form in a temperature range between 273.2 and 272.2 K. In this range, the amount of unfrozen water in a fine-grained soil is still large enough to give a non-negligible permeability. In general, the permeability, K, of a (partially) frozen soil is defined from the water flux equation, Eq. (2), by: In theory, the permeability can be determined by applying a pressure difference on the sample at subzero temperature, T, and measuring the mass flux of water. The sample thickness, Dx, and the pressure difference, Dp, are finite, and the ratio of the two can be regarded as an average value over the state of the sample. Transport is caused by a high pressure at the inlet of the sample, and lower pressure at the sample outlet. The steady water flux should follow from a pressure difference alone. According to Eq. (4), the sample is partially frozen, but isothermal. The experiment requires ice-lens-free samples, otherwise the lens would block the flow. Such an experimental setup was proposed by Black and Miller [2], see the sketch in Fig. 2. Theoretically, it is also possible to follow a different route to access the permeability coefficient, which is provided by Eq. (12) The experiments should be repeated at different temperatures to give K ¼ KðTÞ. Depending on the soil type, even at subzero temperatures, some water can stay in its liquid state in the capillaries (for instance in clay, 30-40% of water content can stay in liquid state at -1°C). In this manner, one can find K ¼ KðTÞ for a possible range of temperatures applicable to a frozen fringe, where unfrozen water and ice are simultaneously present in the pores. The value of K as a function of temperature can be transformed into a function of unfrozen water saturation by means of the freezing characteristic curve. Equation (4) gives Darcy's law for the frozen fringe; and an expression for the first Onsager coefficient:

The thermodynamic frost heave coefficient
We now define the thermodynamic frost heave coefficient that follows from Eq. (2): where T is the average temperature of the frozen fringe. The coefficient s (in kg m À1 K À1 s À1 ) is defined as the ratio of the water flux to the temperature gradient across the frozen fringe at uniform liquid pressure, p. The coefficient s is related to the coupling coefficient L wq in the original set of equations. For a direct measurement of s, we need to measure the temperatures at the boundaries of the frozen fringe and the water flux, at a constant pressure. It is not easy to control the pressure in a partially frozen zone like frozen fringe. An experiment to measure the ratio in Eq. (7) without pore pressure control was called the frost heave test; see e.g. [18], and the coefficient was called the segregation potential SP. We propose here to replace the determination of SP, with a more precise coefficient s. A frost heave cell can be used to determine the thermodynamic frost heave coefficient (Fig. 3). In this experiment, a pretreated, consolidated, homogeneous sample is thermostatted at both ends, and the temperature profile along the sample is recorded. We apply a temperature gradient over the sample length, Dl. The sample has free access to water from the warm end. At steady state, we measure the length of the frozen fringe, Dx, the temperature difference over the frozen fringe and the flux of mass into the system. Ideally, the liquid pressure should be constant in this experiment. Then, s is well-defined as a function of state variables like temperature, while SP depends on pressure gradient in addition to the state variables.
Theory allows us to, alternatively, determine the coefficient s from the permeability test described in Sect. 3.1. By dividing Eqs. (1) and (2) at an isothermal condition, we obtain: The ratio is called the heat of transfer. The energy balance across the frozen fringe provides the insight into the order of magnitude of this quantity. The heat flux out of the frozen layer minus the heat flux into the frozen fringe is the product of the mass flux, _ w, times the enthalpy of freezing, D w;i H(in J kg À1 ). Therefore, the heat of transfer q* must be of the same order of magnitude as D w;i H. It was shown by Kjelstrup and Bedeaux [14] that the heat of transfer is a fraction of the enthalpy difference: where D w;i H ¼ H i À H w , H w is the specific enthalpy of water and H i is the specific enthalpy of ice, and j is a dimensionless number between 0 and 1. The factor depends on the ratio of the soil thermal resistivity to the ice interface resistivity [14]. For instance, in Devon silt, the thermal conductivity k (the inverse of the resistivity) varies between 1.76 (frozen silt) and 1.55 (partially frozen silt) W m À1 K À1 [19], giving a mean value of 1.6 W m À1 K À1 . The conductivity of ice at 0°C is about 2.2 W m À1 K À1 . Using the formula of Kjelstrup and Bedeaux [14], we obtain: The meaning of this number is that 60% of the heat liberated by the phase transition is delivered to the ice, the remaining part is delivered to the soil matrix. Derjaguin and Churaev [7] allocated the whole enthalpy of freezing to the heat of transfer, meaning that all heat was delivered to the ice. By introducing Eq. (9) into Eq. (8) we obtain the coupling coefficient: The thermodynamic frost heave coefficient, s, can now be related to permeability by: The estimate of j can be evaluated by measuring s, and permeability K. Note that the permeability, K, refers to the temperature, T, and so will s. The other terms are constants. It is possible to perform the measurement defined by Eq. (7) than by Eq. (4). Then, Eq. (12) provides an alternative route to permeability.

Thermal conductivity
The thermal conductivity, k, is defined from the heat flux Eq. (1) at uniform pressure: The thermal conductivity (in W m À1 K À1 ) is also a function of unfrozen water saturation. From Eq. (13), we can determine the thermal conductivity, when the pore water pressure is constant along the frozen fringe.
Here, we can use the experimental setup proposed by Côté and Konrad [4]. In this experiment, the pretreated, consolidated sample is mounted between two temperaturecontrolled plates. Thermo-electrical heat flux meter (flux plates) equipped with thermocouples is placed between the sample and the plates. The temperatures for the top and bottom plates are set using thermal baths; and measured using thermocouples. The laboratory setup is illustrated in Fig. 4. As unfrozen water can play an important role in frozen states, thermal conductivities at different freezing gradients can be measured to estimate the thermal conductivity as a function of temperature or unfrozen water saturation.

The pressure difference during water transport across the fringe
Knowing L ww and L wq (as described in Sects. 3.1 and 3.2), we can calculate the pressure difference across the fringe. This expression is relevant when there is a net flux of water to the ice (a growing ice lens):

A practical set of flux-force relations
We can now contract all the information discussed above, introduce the measured properties, and obtain a practical set of flux equations: In the above equations, we have assumed T=T 0 % 1. The last equation is interesting because it shows that a maximum pressure difference exists for a given temperature difference. A temperature difference across a frost heave susceptible frozen soil leads to a water flux. According to Eq. (14), a counter pressure will build, until in the end the thermal force is balanced by the pressure force. At this point, when the pressure difference reaches its maximum value, and the balance of forces is achieved, the water flux will stop. The state, described by the balance of forces, is called the Soret equilibrium state. The ratio of the driving forces in Soret equilibrium becomes: The balance of forces is illustrated in Fig. 5, for j ¼ 1, DT ¼ 0:5°C and s=K from Eq. (12). Figure shows a linear reduction in the temperature from the top boundary of the unfrozen layer to the bottom of the ice lens. Counter to that is a gradual linear increase in the water pressure, ending with a maximum value at the ice. In this situation, there is no net driving force for water flux in Eq. (16).
The maximum frost heave pressure is given by the thermodynamic frost heave coefficient and the permeability. One may expect that these properties depend on the soil in question. Takashi et al. [26] reported a limiting value of 1.14 MPa/K for Manaita bridge clay and Negishei silt, which is corresponding to j ¼ 1 in the present theory. The value was considerably reduced as the temperature of the cooling plate went down. Interestingly, the authors reported that it was possible to move from one Soret equilibrium to another, following the same slope. This is expected for a reversible phenomenon and supports the above description.
Equations (15) and (16), together with equations for mass and energy conservation (given by Eqs. 18 and 20 in ''Appendix''), were used to model the final stage of a frost heave test on Devon silt, using data reported by Konrad and Morgenstern [18]. The onset of formation of the final ice lens was considered as our time reference (t ¼ 0) in this simulation, and the amount of heave from this point was calculated and compared with the experimental result. In this test, the soil sample was initially unfrozen, with a temperature of ? 3°C. The top surface temperature was next reduced to -5.5°C to trigger the freezing process. It remained constant during the test. Water was freely available at the base at atmospheric pressure. The initial porosity of the 78 mm soil sample was 38%. The hydraulic conductivity of the unfrozen soil was 1 9 10 -7 cm/s, the soil freezing characteristic curve and relative permeability relation with temperature are available in Azmatch [1]. The thermal conductivity of unfrozen, partially frozen and fully frozen states were 1.47, 1.51 and 1.76 W/m o C, respectively. The thermodynamic frost heave coefficient for this analysis is calculated from Eq. (12) as s=K ¼ 1:22 Â 10 9 j (kg 2 m À3 K À1 s À2 ). Frost heave is computed for different values of j ¼ 1, 0:9, 0:6 and presented in Fig. 6. As shown in figure, the result corresponding to j ¼ 0:9 has the best match with the experiment.
The water pressure and temperature profiles at the end of the experiment were computed from the model and are shown in Fig. 7. Below the freezing front, a suction of  Comparison between the measured and predicted frost heave for Devon silt around 27 kPa is predicted by our model. This suction value is in the range of those reported on Devon silt by Konrad and Duquennoi [16]. Note that at this stage, temporal variation of the temperature profile, and thus the pressure, are extremely slow. We might practically say that the system has reached a steady state condition. However, this is not completely true, since the growing ice lens will increase the length of the sample and the distance of the frozen fringe from the cooling plate. Thus, the thermal condition in the frozen fringe gets slowly affected.
The model was also used to simulate the steady state part of the frost heave test on Xuzhou silty clay conducted by Zhou and Zhou [32]. The 13 cm soil sample in this test was initially unfrozen, with a temperature of ? 6°C. The bottom end of the column was maintained unchanged at this temperature, while the top end was cooled and kept at -12°C. The sample was allowed to take up water from the bottom, at atmospheric pressure. The test was conducted in an open system without any overburden. The initial porosity of the sample was 30%. The permeability of the unfrozen soil was 1 9 10 -8 cm/s, the soil freezing characteristic curve, relative permeability and thermal conductivity relations with temperature are given in Zhou and Zhou [32]. The thermodynamic frost heave coefficient for this analysis is again calculated from Eq. (12) as s=K ¼ 7:32 Â 10 8 (kg 2 m À3 K À1 s À2 ); i.e. j ¼ 0:6. Simulation results are compared with the test data in Fig. 8, and reasonable agreement has been achieved.

Conclusion
We have developed a discrete description of the coupled transport of heat and mass that takes place across the frozen fringe when the temperature sinks below the freezing temperature in frost heave susceptible soils. The description can be seen as an extension of the works of Derjaguin and Churaev [6], Førland and Kjelstrup Ratkje [8] and Kjelstrup and Bedeaux [14], and treats the fringe in a discrete manner, with interfaces controlled by the table of water and the growing ice lens. In this domain, Darcy's and Fourier's laws do not apply. An emphasis was made on how to measure the transport coefficients. A new coefficient was defined; the thermodynamic frost heave coefficient, expressing the ability to increase the pressure at the ice lens. The coefficient is similar, yet different from the segregation potential defined by Konrad and Morgenstern [18]. The transport coefficients were defined in Eqs. (4), (7) and (13), and the Onsager symmetry relation was applied to find the expression for the maximum frost heave pressure, characteristic for each soil. The analysis was illustrated with numerical examples from Devon silt and Xuzhou silty clay. This completes the thermodynamic analysis of the transport phenomena in frost heave. A mechanical analysis will follow in upcoming works.

Availability of data and material
All the used data were cited, and the relevant references can be found in the reference list.

Appendix: equations of mass and energy conservation
The mass conservation equation for water and ice phases can be written as: where / denotes the porosity of the soil. The following boundary condition is applied to the mass conservation equation at the ice lens (interface 2-3 in Fig. 1): where _ h is the heave rate. Assuming local thermal equilibrium, a single energy balance equation can be used to describe the heat transfer process in a multiphase system. This assumption implies that all the phases at each spatial point reach thermal equilibrium instantaneously together. Neglecting the kinetic energy, viscous and intrinsic dissipation, the energy conservation equation can be written as: where C is the effective heat capacity of frozen soil, V denotes the specific volume of frozen soil, C w is the heat capacity of water and _ m w!i is standing for the phase change of water to ice. Code availability The in-house code used for solving the equations is available upon request.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.