1 Introduction

Approximately 17% of the northern hemisphere is permanently frozen ground (permafrost) and another 34% are seasonally frozen soils [42]. In China, seasonal frozen soils occupy an area of 513.7 km2, which is 53.5% of the total area of the 13 related provinces’ municipalities and autonomous regions. The areas underlain by permafrost and seasonally frozen soils experience serious engineering problems. Statistics indicate that 83% of the canals in the highly irrigated district of Heilongjiang have suffered frost damage due to inappropriate anti-frost heaving measures. In Jilin and northern Xinjiang, more than half of the concrete trunk and branch canals exhibit different degrees of frost heave damage. In the Hetao region of Inner Mongolia, approximately 70% of the hydraulic buildings exhibit various degrees of freezing and thawing damage. In the Shule area of the Gansu Province, before repair operations, 54% of the canals had serious frost heave damage and lost their anti-seepage properties. Therefore, safe operation of the canals cannot be guaranteed. Both construction and safe operations of these systems in cold areas are severely restricted by permafrost foundations. For example, before 2015, frost-induced uplifting of foundation piles resulted in the Harbin–Dalian and Panjin–Yingkou high-speed railways reducing their operating speeds to 200 km/h in winter, which is considerably lower than the standard speed of 300–350 km/h. The physical and mechanical properties of frozen soil and frost heave models must be thoroughly understood for construction projects in cold regions.

There are two main reasons that cause damage of structures in seasonally frozen ground: the heaving and thawing settlement of the foundation soil. Owing to wider capillary water migration under the same temperature and moisture conditions, fine-grained soil has a higher frost heave sensitivity than other soil types. In addition, a recent study [25] has demonstrated that with a sufficient water supply, significant frost heave can also occur in gravel with a limited amount of fine material. Frost heaving and freezing in soils are complex physical processes. Frost damage appears during repeated freezing and thawing cycles of the soil throughout the engineering life cycle. In particular, frost damage is the result of long-term repeated interactions between the soil and a superstructure. In this process, soil heaving is caused by phase transitions of the water in the soil from fluid to solid, but settlement is due to reverse phase transition. For an accurate analysis of the interaction between permafrost and buildings, accurate models of permafrost constitutive and frost heave are required. A frost heave model considers various factors, such as soil hydrodynamics, the physics of frozen soils, frozen soil mechanics, thermodynamics, and engineering mechanics. For example, the rigid ice model [28, 32] and the segregation potential model [12, 13, 27] can accurately explain the formation of an ice lens. However, these models are primarily used to simulate one-dimensional freezing of saturated homogeneous soils. The dissipation of an ice lens in the rigid ice model has been rarely reported. In addition, the interaction between frozen soils and building has rarely been considered. The HFHM [5, 6, 8, 21, 23, 34], also known as the thermal–moisture–deformation coupled model or physical field model, cannot adequately describe the characteristics of the generation and elimination of the ice lens, but they have an advantage in the determination of the macro-characteristics of the temperature distribution, water migration, frost heave, and thaw settlement. Numerous indoor and field monitoring results [8, 9, 11, 16, 31] have shown that HFHM is suitable to describe the energy and mass transfers during the freezing and thawing processes. In addition, the physical field model also can be used to describe the coupled process in the frozen fringe. A separate-ice frost heave model, proposed by Zhou et al. [44], can describe the growth of a single ice lens by extending the thermal–moisture model and the formation of a new ice lens using Gilpin’s theory. Furthermore, HFHM is highly compatible both with the nonlinear dynamics model [15, 40] and with the contact model [20, 36, 43] for permafrost. A centrifuge model and numerical modeling of a cold-region canal were conducted by Li et al. [19]. The models indicated that the HFHM-computed temperature perfectly agreed with the measured results, but had less accurate results regarding deformation along the lining. This was likely due to complicated multidirectional freezing.

Generally, the HFHM assumes that frost heave occurs when the total water content (ice content plus unfrozen water content) of the feature unit is larger than the critical pore volume [5, 10, 18]. Deformation of the unit is then assumed to be evenly distributed in all of the dimensions and independent of the temperature gradient. The constitutive and frost heave models are considered to be ISO models. The empirical formula for the elastic constants of the frozen soils can be obtained according to the composition of the physical components using the mixed law [24]. The aforementioned two assumptions imply that the models ignore the directionality and distribution of the growth of the ice lens under a temperature gradient, which conflicts with the freezing and frosting evolution and microstructure of the frozen soil. Therefore, a realistic coupled thermal–moisture–deformation model is required to accurately reflect the composition, mesoscopic structure, and frost heave evolution. Additionally, the model should be able to describe the orthogonality of the frost heave and the microstructure of the frozen soil in the direction of the temperature gradient and frozen fronts. An appropriate coupled thermal–moisture–deformation model needs to contain a constitutive model and a model for frost heave in frozen soils.

The ISO constitutive model has been widely used. However, it is inconsistent with the profile of layered ice lenses. The orthotropic constitutive model is more realistic than the ISO constitutive model. However, it is complicated by having nine elastic constants, which cannot be easily determined using tests. Transversely isotropic (TISO) constitutive models combine the advantages of the aforementioned two models. TISO models have been widely used. Studies have indicated that fine-grained frozen soils belong to TISO materials. Most frozen soil tests involve triaxial tests with a constant temperature and overburden without water supply, creep tests, and frost heave tests. Up to now, there has been no feasible method or device that can directly determine the TISO elastic constants for frozen soils. Although a few mechanical parameters can be obtained using tests, the testing process is too complicated to effectively utilize in multiphysical coupled simulations. Therefore, a TISO constitutive model has been proposed for the theoretical calculation of the elastic constants of frozen soils [38]. However, the constitutive model does not consider the static balance between the ice and soil and does not take into account their deformation compatibility. Furthermore, the TISO constitutive models are seldom built for frozen soils. Even though plastic or creep models for ISO frozen soils have been studied since the 1970s [14], related models for TISO frozen soils have not been reported.

In this study, the serial model and parallel model are proposed on the basis of both the mesoscopic profile of fine-grained frozen soils and the theory of composite mechanics. This is done by assuming an ideal bond between the ice layers and soil layers. The elastic constants of TISO frozen soils with ice lenses are derived by combining the primary equations of Hooke’s law, equilibrium equations, and strain compatibility equations. In addition, in terms of the strength index of frozen soil at the principle directions, an elastoplastic model for TISO frozen soils is proposed based on the Hill orthotropic plasticity theory. Furthermore, a modified HFHM for freezing soils is presented that considers the evolution of the thermal–moisture–deformation interaction and the TISO characteristics of freezing soil.

2 Frost heave characteristics of freezing fine-grained soils

The results of the indoor unidirectional freezing test [33] indicated that when the freezing front is irregularly developing during freezing, the ice lenses and the dehydrated soil layers form laminated structures in the frozen zone on the mesoscale (Fig. 1). Furthermore, soil particles, unfrozen water, crystal ice, air, and small serial cracks were present in the frozen zone because of the complex formation process of the fine-grained soils. Therefore, obvious macroscopic differences were observed in the linear and nonlinear mechanical performances of the frozen soils along their principal directions. Arenson et al. [1] considered that fine-grained soils had a clear stratification in the frozen zone when the ice content was between 10 and 90%. Results from a series of unconfined unidirectional freezing tests and an expanded foundation model test [37] indicated that the frozen silt could be modeled as TISO under natural freezing conditions. Correspondingly, the nine elastic constants of the orthotropic constitutive model can be reduced to five in the TISO constitutive model. Moreover, Wang et al. [39] developed an empirically nonlinear constitutive model for TISO frozen silt on the basis of the axial and lateral deformations observed in unconfined uniaxial compression tests, which were performed at temperatures ranging from −9 to −5 °C.

Fig. 1
figure 1

Snapshots of a cylinder sample of frozen clay in a freezing test by Taber [33], indicating that the ice lenses and dehydrated soil layers formed alternately laminated structures during freezing. a Clay frozen without overburden. b Clay frozen under pressure of 15 kilograms per square centimeter

Full size image

The physical properties of a frozen soil depend on the freezing process. Uniaxial freezing tests of frost-susceptible soils indicate that migration of the moisture in fine-grained soils exhibits directional dependence during freezing. Initially, due to a high temperature gradient near the cold boundary, the frozen fringe advances quickly from the cool side to the warm side. Rapid cooling causes in situ freezing of the water in soil. Heaving results in water migration being impeded by a dispersive ice crystal. As the temperature gradually decreases, the frozen front acts as a sink. Water migrates continuously if a temperature gradient exists and an external water supply is available. If the sum of the unfrozen water and accumulated ice exceeds the pore volume in the soil, layers of clear ice may form parallel to the surface due to considerable frost heaving. Moreover, unconfined uniaxial frost heave tests of fine-grained soil indicate that thermal expansion along the temperature gradient direction is larger than that along the perpendicular direction. Thus, the frozen fine-grained soil deforms in a TISO manner rather than in an ISO manner, which determines different macroscopically mechanical properties and frost heave. Therefore, the TISO constitutive model is more reasonable than the ISO constitutive model for describing both frost heaving and the interaction between frozen soil and a foundation.

To predict the mechanical behavior of a TISO frozen soil in a construction project, the material constants must be first obtained. However, few laboratory tests have been conducted on the mechanics of an ice-embedded frozen soil due to such reasons as inhomogeneous distribution, environmental sensitivity, and large uncertainties of the basic material properties. An analytical method for determining the elastic constants of frozen soil with a lenticular structure is presented in this paper.

3 Transversely isotropic elastic constants for frozen soils with ice lenses

The elastic moduli and Poisson’s ratios can be used to determine the characteristic features of the freezing processes of frost-susceptible soils. Ice provides mechanical bonds on the interface between a soil particle and the bedding planes. The ice lenses and soil layers in the representative volume element (RVE) are assumed to possess deformability. Furthermore, the soil with lenticular structure can be considered as a single-layered composite when the REV is suitably large. In the local Cartesian coordinate system, anisotropy of material occurs in the 2-direction, that is, the direction of the heat flow lines (Fig. 2). The principal material axes are aligned with the coordinate axes, respectively. Material in planes 1–3 shows isotropic (ISO) behavior. The strain vector is defined as ε = [ε11, ε22, ε33, γ12, γ23, γ31]T, and the stress vector is defined as σ = [σ1, σ2, σ3, σ11, τ12, τ23, τ31]T. εii and σi represent the normal strain and stress along the i-direction, respectively. γij and τij represent the shear strain and stress, respectively.

Fig. 2
figure 2

RVE schematic of a single-layered composite model for frozen soil. a Serial model for estimating elastic constants E2 and υ21 (= υ23). b Parallel model for estimating elastic constants E1, υ12, and υ13. RVE representative volume element

Full size image

3.1 Anisotropic character of frost heave

The frost heave of frozen fine-grained soil exhibited temperature dependence due to the formation of the ice lens [22, 37]. Hence, a dimensionless quantity, ξ (the frost heave partition coefficient), is used to illustrate the anisotropy of the volumetric growth. The plastic strain increment can be defined as follows:

$${\varvec{\upvarepsilon}}_{{{\text{vh}},123}} = \left( {\theta_{\text{w}} + \theta_{\text{i}} - n_{0} } \right)\left[ {\begin{array}{*{20}c} {\left( {1 - \xi } \right)/2} & 0 & 0 \\ 0 & \xi & 0 \\ 0 & 0 & {\left( {1 - \xi } \right)/2} \\ \end{array} } \right],$$
(1)

where εvh,123 is the plastic strain tensor for the entire frozen composite due to volumetric growth in the local system. The volumetric contents of unfrozen water and ice and the initial porosity are represented as θw, θi, and θ0, respectively. The value of ξ ranges from 0.33 to 1.0 for the ISO growth and unidirectional growth of ice lenses.

3.2 Stiffness matrix of frozen soil

Conventionally, the five independent elastic constants in the TISO constitutive model are the Young’s modulus and Poisson’s ratio in planes 1–3 (E1 and υ13, respectively), Young’s modulus and Poisson’s ratio in the 2-direction (E2 and υ12, respectively), and the shear modulus in the 2-direction (G12). The stiffness matrix reads:

$${\mathbf{D}}_{123} = \left[ {\begin{array}{*{20}l} {d_{11} } \hfill & {d_{12} } \hfill & {d_{13} } \hfill & {} \hfill & {} \hfill & {} \hfill \\ {} \hfill & {d_{22} } \hfill & {d_{12} } \hfill & {} \hfill & 0 \hfill & {} \hfill \\ {} \hfill & {} \hfill & {d_{11} } \hfill & {} \hfill & {} \hfill & {} \hfill \\ {} \hfill & {} \hfill & {} \hfill & {d_{44} } \hfill & {} \hfill & {} \hfill \\ {} \hfill & {{\mathbf{sym}}} \hfill & {} \hfill & {} \hfill & {\left( {d_{11} - d_{13} } \right)/2} \hfill & {} \hfill \\ {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {d_{44} } \hfill \\ \end{array} } \right],$$
(2)

where

$$\begin{aligned} & d_{11} = \frac{{1 - \upsilon_{12} \upsilon_{21} }}{{\left( {1 - \upsilon_{13} - 2\upsilon_{21} \upsilon_{12} } \right)\left( {1 + \upsilon_{13} } \right)}}E_{1} , \\ & d_{22} = \frac{{1 - \upsilon_{13} }}{{1 - \upsilon_{13} - 2\upsilon_{21} \upsilon_{12} }}E_{2} , \\ & d_{12} = \frac{{\upsilon_{21} + \upsilon_{13} \upsilon_{21} }}{{\left( {1 - \upsilon_{13} - 2\upsilon_{21} \upsilon_{12} } \right)\left( {1 + \upsilon_{13} } \right)}}E_{1} , \\ & d_{13} = \frac{{\upsilon_{13} + \upsilon_{21} \upsilon_{12} }}{{\left( {1 - \upsilon_{13} - 2\upsilon_{21} \upsilon_{12} } \right)\left( {1 + \upsilon_{13} } \right)}}E_{1} , \\ & d_{44} = G_{12} , \\ & \upsilon_{12} /E_{1} = \upsilon_{21} /E_{2} . \\ \end{aligned}$$

By assuming deformability, the serial model and the parallel model [45] can be used to estimate the elastic constants and unbalanced shear force between the ice and soil layers due to different Poisson’s ratios of contents. The elastic modulus and Poisson’s ratio for the ice and soil matrix are represented as Ei, υi, Em, and υm, respectively. The volume of the ice lenses and soil in an RVE are represented as θi and θm = 1 − θi, respectively. The stresses at the interface are represented as σji and σjm (j = 1, 2, 3 are the principal material axes in Fig. 2).

3.2.1 Serial model

A modified serial model that considers the deformability and the interlayer interaction is used to estimate the elastic constants of the frozen soil along the 2-direction (Fig. 2a). The RVE is assumed to be subjected to an axial load, σ2, along the 2-direction. The equilibrium along the 2-direction suggests σ2i = σ2m. The frictional force between the layers is absent if the deformation of an individual layer is either independent or the layers have the same Young’s modulus and Poisson’s ratio. Then, the effective elastic constants of the RVE can be obtained using the general rule of mixtures. For consistent lateral deformation of the two layers, the interactions parallel to the interface (σ1i, σ3i, σ1m, and σ3m) are nonzero because of the interfacial adhesion and the differences between the Poisson’s ratios for ice and soil. Considering symmetry in the planes 1–3, then:

$$\begin{aligned} \sigma_{{1{\text{m}}}} & = \sigma_{{3{\text{m}}}} = \sigma_{\text{m}} \\ \sigma_{{1{\text{i}}}} & = \sigma_{{3{\text{i}}}} = \sigma_{\text{i}} . \\ \end{aligned}$$
(3)

The equilibrium equations and deformability along the principal axes yield:

$$\begin{aligned} & \sigma_{{2{\text{i}}}} = \sigma_{{ 2 {\text{m}}}} = \sigma_{2} \\ & \sigma_{\text{m}} \left( {1 - \theta_{\text{i}} } \right) + \sigma_{\text{i}} \theta_{\text{i}} = 0, \\ \end{aligned}$$
(4)
$$\begin{aligned} \varepsilon_{1} & = \varepsilon_{3} = \frac{1}{{E_{\text{i}} }}\left( {\sigma_{\text{i}} - \nu_{\text{i}} \sigma_{\text{i}} - \nu_{\text{i}} \sigma_{2} } \right) \\ & = \frac{1}{{E_{\text{m}} }}\left( {\sigma_{\text{m}} - \nu_{\text{m}} \sigma_{\text{m}} - \nu_{\text{m}} \sigma_{2} } \right). \\ \end{aligned}$$
(5)

Hence,

$$\sigma_{\text{m}} = \frac{{\left( {\nu_{\text{m}} - \nu_{\text{i}} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)\sigma_{2} }}{{1 - \upsilon_{\text{m}} + \left( {1 - \upsilon_{\text{i}} } \right)\frac{{1 - \theta_{\text{i}} }}{{\theta_{i} }}\frac{{E_{\text{m}} }}{{E_{\text{i}} }}}},$$
(6)
$$\sigma_{\text{i}} = - \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }}\sigma_{\text{m}} = - \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }}\frac{{\left( {\nu_{\text{m}} - \nu_{\text{i}} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)\sigma_{2} }}{{1 - \upsilon_{\text{m}} + \left( {1 - \upsilon_{\text{i}} } \right)\frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }}\frac{{E_{\text{m}} }}{{E_{\text{i}} }}}}.$$
(7)

For a cubic element within the calculated domain, Hooke’s law reads:

$$\varepsilon_{2} = \frac{1}{{E_{\text{i}} }}\left( {\sigma_{2} - 2\upsilon_{\text{i}} \sigma_{\text{i}} } \right)\theta_{\text{i}} + \frac{1}{{E_{\text{m}} }}\left( {\sigma_{2} - 2\upsilon_{\text{m}} \sigma_{\text{m}} } \right)\left( {1 - \theta_{\text{i}} } \right).$$
(8)

By substituting Eqs. 6 and 7 into Eqs. 5 and 8, the strain components of the RVE along the principal axes can be represented using the elastic constants of soil and ice. Hence:

$$\varepsilon_{1} = \varepsilon_{3} = \frac{{ - \sigma_{2} }}{{E_{\text{i}} }}\frac{{\left( {1 - \nu_{\text{m}} } \right)\nu_{\text{i}} + \nu_{\text{m}} \left( {1 - \nu_{\text{i}} } \right)\frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }}}}{{\left( {1 - \nu_{\text{m}} } \right) + \left( {1 - \nu_{\text{i}} } \right)\frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}}},$$
(9)
$$\varepsilon_{2} = \frac{{\sigma_{2} }}{{E_{\text{m}} }}\left[ {\left( {1 - \theta_{\text{i}} + \theta_{\text{i}} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right) - \frac{{2\left( {1 - \theta_{\text{i}} } \right)\left( {\nu_{\text{m}} - \nu_{\text{i}} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} }}{{\left( {1 - \nu_{\text{m}} } \right) + \left( {1 - \nu_{\text{i}} } \right)\frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}}}} \right],$$
(10)

and the two elastic constants of the element, E2 and υ21, can be obtained. They are:

$$\begin{array}{*{20}c} {E_{2} = \sigma_{2} /\varepsilon_{2} = E_{\text{m}} \left[ {\left( {1 - \theta_{\text{i}} + \theta_{\text{i}} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right) - \frac{{2\left( {1 - \theta_{\text{i}} } \right)\left( {\nu_{\text{m}} - \nu_{\text{i}} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} }}{{\left( {1 - \nu_{\text{m}} } \right) + \left( {1 - \nu_{\text{i}} } \right)\frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}}}} \right]^{ - 1} } \\ \end{array} ,$$
(11)
$$\begin{array}{l} {\nu _{{21}} = \nu _{{23}} = - \varepsilon _{3} /\varepsilon _{2} } \hfill \\ { = \frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}\left[ {\nu _{{\text{i}}} \left( {1 - \nu _{{\text{m}}} } \right) + \nu _{{\text{m}}} \left( {1 - \nu _{{\text{i}}} } \right)\frac{{1 - \theta _{{\text{i}}} }}{{\theta _{{\text{i}}} }}} \right]} \hfill \\ {\left\{ {\left( {1 - \theta _{{\text{i}}} + \theta _{{\text{i}}} \frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right)} \right.} \hfill \\ {\left[ {\left( {1 - \nu _{{\text{m}}} } \right) + \left( {1 - \nu _{{\text{i}}} } \right)\frac{{1 - \theta _{i} }}{{\theta _{{\text{i}}} }}~\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right]} \hfill \\ {\left. { - 2\left( {1 - \theta _{{\text{i}}} } \right)\left( {\nu _{{\text{m}}} - \nu _{{\text{i}}} \frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right)^{2} } \right\}^{{ - 1}} }.\end{array}$$
(12)

3.2.2 Parallel model

The parallel model accounts for deformability and is suitable for predicting the linear elasticity of frozen soil in the symmetric plane (planes 1–3 in Fig. 2b). The RVE is assumed to be subjected to an axial load, σ1, along the 1-direction. The equilibrium equations and deformability along the principal axes satisfy the following:

$$\begin{aligned} & \sigma_{{1{\text{i}}}} \theta_{\text{i}} + \sigma_{{1{\text{m}}}} \left( {1 - \theta_{\text{i}} } \right) = \sigma_{1} \cdot 1 \\ & \sigma_{{3{\text{i}}}} \theta_{\text{i}} + \sigma_{{3{\text{m}}}} \left( {1 - \theta_{\text{i}} } \right) = 0 \\ \end{aligned}$$
(13)
$$\begin{aligned} \varepsilon_{1} & = \frac{{\sigma_{{1{\text{i}}}} - \nu_{\text{i}} \sigma_{{3{\text{i}}}} }}{{E_{\text{i}} }} = \frac{{\sigma_{{1{\text{m}}}} - \nu_{\text{m}} \sigma_{{3{\text{m}}}} }}{{E_{\text{m}} }} \\ \varepsilon_{3} & = \frac{{\sigma_{{3{\text{i}}}} - \nu_{\text{i}} \sigma_{{1{\text{i}}}} }}{{E_{\text{i}} }} = \frac{{\sigma_{{3{\text{m}}}} - \nu_{\text{m}} \sigma_{{1{\text{m}}}} }}{{E_{\text{m}} }}. \\ \end{aligned}$$
(14)

Hooke’s law can be expressed as:

$$\varepsilon_{2} = - \frac{{\nu_{\text{i}} \theta_{\text{i}} \left( {\sigma_{{1{\text{i}}}} + \sigma_{{3{\text{i}}}} } \right)}}{{E_{\text{i}} }} - \frac{{\nu_{\text{m}} \left( {1 - \theta_{\text{i}} } \right)\left( {\sigma_{{1{\text{m}}}} + \sigma_{{3{\text{m}}}} } \right)}}{{E_{\text{m}} }}.$$
(15)

By combining Eqs. 13–15, we have:

$$\sigma_{{1{\text{i}}}} = \frac{{\sigma_{1} }}{{\theta_{\text{i}} }}\left[ {1 - \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }}\frac{{E_{\text{m}} }}{{E_{\text{i}} }}\frac{{1 + \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }} - \nu_{i} \left( {\nu_{\text{m}} + \nu_{\text{i}} \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)}}{{\left( {1 + \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} - \left( {\nu_{\text{m}} + \nu_{\text{i}} \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} }}} \right],$$
(16)
$$\sigma_{{1{\text{m}}}} = \left( {\sigma_{1} \cdot 1 - \sigma_{{1{\text{i}}}} \theta_{i} } \right)/\left( {1 - \theta_{\text{i}} } \right),$$
(17)
$$\sigma_{{3{\text{i}}}} = - \frac{{\sigma_{1} }}{{\theta_{\text{i}} }}\frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }}\frac{{E_{\text{m}} }}{{E_{\text{i}} }}\frac{{\nu_{\text{m}} - \nu_{\text{i}} }}{{\left( {1 + \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} - \left( {\nu_{\text{m}} + \nu_{\text{i}} \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} }},$$
(18)
$$\sigma_{3m} = - \sigma_{3i} \theta_{i} /\left( {1 - \theta_{i} } \right).$$
(19)

By substituting Eqs. 16–19 into Eqs. 14 and 15, we have:

$$\varepsilon_{1} = \frac{{\sigma_{1} }}{{\theta_{\text{i}} E_{\text{i}} }}\frac{{1 - \nu_{\text{m}}^{2} + \left( {1 - \nu_{\text{i}}^{2} } \right)\frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}}}{{\left( {1 + \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} - \left( {\nu_{\text{m}} + \nu_{\text{i}} \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} }},$$
(20)
$$\varepsilon _{2} = - \frac{{\sigma _{1} }}{{\theta _{{\text{i}}} E_{{\text{i}}} }}\left\{ {\nu _{{\text{i}}} \theta _{{\text{i}}} + \frac{{\left( {1 - \theta _{{\text{i}}} } \right)\left( {\nu _{{\text{m}}} - \nu _{{\text{i}}} \frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right)\left( {1 - \nu _{{\text{i}}} } \right)}}{{\left( {1 + \frac{{1 - \theta _{{\text{i}}} }}{{\theta _{{\text{i}}} }}~\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right)^{2} - \left( {\nu _{{\text{m}}} + \nu _{{\text{i}}} \frac{{1 - \theta _{{\text{i}}} }}{{\theta _{{\text{i}}} }}~\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right)^{2} }}\left. {\frac{{\left[ {1 + \nu _{{\text{m}}} + \left( {1 + \nu _{{\text{i}}} } \right)\frac{{1 - \theta _{{\text{i}}} }}{{\theta _{{\text{i}}} }}~\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right]}}{{\left( {1 + \frac{{1 - \theta _{{\text{i}}} }}{{\theta _{{\text{i}}} }}~\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right)^{2} - \left( {\nu _{{\text{m}}} + \nu _{{\text{i}}} \frac{{1 - \theta _{{\text{i}}} }}{{\theta _{{\text{i}}} }}~\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}} \right)^{2} }}} \right\}} \right.$$
(21)
$$\varepsilon_{3} = \frac{{\sigma_{1} }}{{\theta_{\text{i}} E_{\text{i}} }}\frac{{ - \nu_{\text{m}} \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }} - \nu_{\text{i}} \left( {1 - \nu_{\text{m}} - \nu_{\text{i}} \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)}}{{\left( {1 + \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} - \left( {\nu_{\text{m}} + \nu_{\text{i}} \frac{{1 - \theta_{\text{i}} }}{{\theta_{\text{i}} }} \frac{{E_{\text{m}} }}{{E_{\text{i}} }}} \right)^{2} }}.$$
(22)

The three elastic constants, E1, υ12, and υ13, can be also obtained.

$$\begin{aligned} E_{1} & = \sigma _{1} /\varepsilon _{1} \\ & = \left[ {E_{{\text{i}}} \theta _{{\text{i}}} + E_{{\text{m}}} \left( {1 - \theta _{{\text{i}}} } \right)} \right] \\ & \quad + \frac{{\left( {\nu _{{\text{m}}} - \nu _{{\text{i}}} } \right)^{2} \theta _{{\text{i}}} \left( {1 - \theta _{{\text{i}}} } \right)E_{{\text{i}}} E_{{\text{m}}} }}{{\left( {1 - \nu _{{\text{m}}}^{2} } \right)\theta _{{\text{i}}} E_{{\text{i}}} + \left( {1 - \nu _{{\text{i}}}^{2} } \right)\left( {1 - \theta _{{\text{i}}} } \right)E_{{\text{m}}} }}, \\ \end{aligned}$$
(23)
$$\begin{aligned} \nu _{{12}} & = - \varepsilon _{2} /\varepsilon _{1} \\ & = \left[ {\nu _{{\text{i}}} \theta _{{\text{i}}} + \nu _{{\text{m}}} \left( {1 - \theta _{{\text{i}}} } \right)} \right] + \left( {\nu _{{\text{m}}} - \nu _{{\text{i}}} } \right)\theta _{{\text{i}}} \left( {1 - \theta _{{\text{i}}} } \right) \\ & \quad \frac{{\nu _{{\text{m}}} \left( {1 + \nu _{{\text{m}}} } \right) - \nu _{{\text{i}}} \left( {1 + \nu _{{\text{i}}} } \right)\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}}}{{\left( {1 - \nu _{{\text{m}}}^{2} } \right)\theta _{{\text{i}}} + \left( {1 - \nu _{{\text{i}}}^{2} } \right)\left( {1 - \theta _{{\text{i}}} } \right)\frac{{E_{{\text{m}}} }}{{E_{{\text{i}}} }}}}, \\ \end{aligned}$$
(24)
$$\begin{aligned} \nu_{13} & = - \varepsilon_{3} /\varepsilon_{1} \\ & = \nu_{\rm i} + \frac{{\left( {\nu_{\rm m} - \nu_{\rm i} } \right)\left( {1 - \nu_{\rm i}^{2} } \right)\left( {1 - \theta_{\rm i} } \right)\frac{{E_{\rm m} }}{{E_{\rm i} }}}}{{\left( {1 - \nu_{\rm m}^{2} } \right)\theta_{\rm i} + \left( {1 - \nu_{\rm i}^{2} } \right)\left( {1 - \theta_{\rm i} } \right)\frac{{E_{\rm m} }}{{E_{\rm i} }}}}. \\ \end{aligned}$$
(25)

3.2.3 Shear modulus

It is assumed that the layers do not separate under shear deformation if an ideal bond exists between them. Thus, the two layers have the same shear strain. Therefore, the equilibrium equation and Hooke’s law can be written as:

$$\tau_{12} = \tau_{\text{i}} \theta_{\text{i}} + \tau_{\text{m}} \theta_{\text{m}} ,$$
(26)
$$\tau_{12} = G_{12} \gamma_{12} ,\quad \tau_{\text{i}} = G_{\text{i}} \gamma_{\text{i}} ,\quad \tau_{\text{m}} = G_{\text{m}} \gamma_{\text{m}} ,$$
(27)

where τ12, γ12, and G12 are the shear stress, engineering shear strain, and effective shear modulus, respectively, in the planes 1–2. Substituting γi = γm = γ12 into Eqs. 26 and 27, G12 reads:

$$G_{12} = G_{\text{i}} \theta_{\text{i}} + G_{\text{m}} \theta_{\text{m}} .$$
(28)

In most analyses of frozen ground engineering, the plane strain model is adopted. The computation will be performed in an arbitrary global coordinate system x–y. The strain increments due to frost heave and the stiffness matrix can be transformed to the x, y coordinate system using the following transformation rule:

$$\begin{aligned} &{\varvec{\upvarepsilon}}_{0} = \left[ {\varepsilon_{x} ,\varepsilon_{y} ,\varepsilon_{z} ,\gamma_{xy} ,\gamma_{yz} ,\gamma_{zx} } \right]^{\text{T}} = {\mathbf{T}}^{ - 1} \left[ {\varepsilon_{0,11} ,\varepsilon_{0,22} ,\varepsilon_{0,33} ,0,0,0 } \right]^{\text{T}} \\ &{\mathbf{D}} = {\mathbf{T}}^{ - 1} {\mathbf{D}}_{123} {\mathbf{T}} \\ &{\mathbf{T}} = \left[ {\begin{array}{*{20}c} {\cos ^{2} \beta } & {\sin ^{2} \beta } & 0 \\ {\sin ^{2} \beta } & {\cos ^{2} \beta } & 0 \\ 0 & 0 & 1 \\ { - 2\cos \beta \sin \beta } & {2\cos \beta \sin \beta } & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} } \right. \\ & \quad \left. {\begin{array}{*{20}c} {\cos \beta \sin \beta } & 0 & 0 \\ { - \cos \beta \sin \beta } & 0 & 0 \\ 0 & 0 & 0 \\ {\cos ^{2} \beta - \sin ^{2} \beta } & 0 & 0 \\ 0 & {\cos \beta } & { - \sin \beta } \\ 0 & {\sin \beta } & {\cos \beta } \\ \end{array} } \right] \\ \end{aligned}$$
(29)

where \(\beta\) is the angle between axis y and the local 2-direction (heat flow direction). It will be synchronously updated as the temperature field changes.

3.3 Ideal plastic model for transversely isotropic frozen soil

In this study, the Hill orthotropic plasticity was applied for describing the elastoplastic behavior of transversely ISO frozen soil. In a complex stress state, the plastic strain rate is generally expressed as [29]:

$${{\dot{\varvec\varepsilon} }_{\bf{p}} }= \lambda \frac{{\partial Q_{\text{p}} }}{{\partial {\varvec{\upsigma}}}}.$$
(30)

A quadratic yield function, Fy (and associated plastic potential, Fy = Qp), in a local coordinate system is given by the orthogonal principal axes [7]:

$$Q_{\text{p}} = F\left( {\sigma_{22} - \sigma_{33} } \right)^{2} + G\left( {\sigma_{33} - \sigma_{11} } \right)^{2} + H\left( {\sigma_{11} - \sigma_{22} } \right)^{2} + 2L\tau_{23}^{2} + 2M\tau_{31}^{2} + 2N\tau_{12}^{2} - 1.$$
(31)

The six parameters F, G, H, L, M, and N are related to the state of anisotropy and can be expressed as:

$$\begin{aligned} F & = \frac{1}{2}\left( {\frac{1}{{\sigma_{{{\text{ys}}2}}^{2} }} + \frac{1}{{\sigma_{{{\text{ys}}3}}^{2} }} - \frac{1}{{\sigma_{{{\text{ys}}1}}^{2} }}} \right) \\ G & = \frac{1}{2}\left( {\frac{1}{{\sigma_{{{\text{ys}}3}}^{2} }} + \frac{1}{{\sigma_{{{\text{ys}}1}}^{2} }} - \frac{1}{{\sigma_{{{\text{ys}}2}}^{2} }}} \right) \\ H & = \frac{1}{2}\left( {\frac{1}{{\sigma_{{{\text{ys}}1}}^{2} }} + \frac{1}{{\sigma_{{{\text{ys}}2}}^{2} }} - \frac{1}{{\sigma_{{{\text{ys}}3}}^{2} }}} \right) \\ L & = \frac{1}{{2\tau_{{{\text{ys}}23}}^{2} }},\quad M = \frac{1}{{2\tau_{{{\text{ys}}31}}^{2} }},\quad N = \frac{1}{{2\tau_{{{\text{ys}}12}}^{2} }}. \\ \end{aligned}$$
(32)

Here, σysi represents the tensile yield stress in the principle direction and τysij represents the yield shear stress on the planes i–j.

For transversely ISO material, F = H and L = N. Hence, the yield function can be rewritten as:

$$\begin{aligned} Q_{\text{p}} & = F\left[ {\left( {\sigma_{22} - \sigma_{33} } \right)^{2} + \left( {\sigma_{11} - \sigma_{22} } \right)^{2} } \right] + G\left( {\sigma_{33} - \sigma_{11} } \right)^{2} + 2L\left( {\sigma_{23}^{2} + \sigma_{12}^{2} } \right) + 2M\sigma_{31}^{2} - 1, \\ F & = \frac{1}{2}\frac{1}{{\sigma_{{{\text{ys}}2}}^{2} }},\quad G = \frac{1}{{\sigma_{{{\text{ys}}1}}^{2} }} - \frac{1}{2}\frac{1}{{\sigma_{{{\text{ys}}2}}^{2} }},\quad L = \frac{1}{{2\sigma_{{{\text{ys}}12}}^{2} }},\quad M = \frac{1}{{2\sigma_{{{\text{ys}}31}}^{2} }}. \\ \end{aligned}$$
(33)

For ISO plasticity, the elastic region satisfies Qp < 0, and the yield surface is defined by Qp = 0.

4 Frost heave modeling with the coupled heat–moisture–deformation model

4.1 Thermal field

The principles of energy, mass, and energy dissipation can be formalized within a multiphase media. Similarly, considering water migration and phase change, the energy conservation principle can be written as:

$$C_{\text{v}} \frac{\partial T}{\partial t} = \nabla \cdot \left( {\lambda \nabla T} \right) + L_{\text{f}} \rho_{\text{i}} \frac{{\partial \theta_{\text{i}} }}{\partial t}.$$
(34)

By neglecting the volumetric heat capacity of air, the semiempirical formulas are used to estimate Cv and λ [4]:

$$\begin{aligned} C_{\text{v}} & = C_{\text{p}} \theta_{\text{p}} + C_{\text{w}} \theta_{\text{w}} + C_{\text{i}} \theta_{\text{i}} \\ \lambda & = \lambda_{\text{p}} \theta_{\text{p}} + \lambda_{\text{w}} \theta_{\text{w}} + \lambda_{\text{i}} \theta_{\text{i}} + \lambda_{\text{a}} \theta_{\text{a}} , \\ \end{aligned}$$
(35)

where T is the temperature; Cv is the effective volumetric heat capacity; λ is the thermal conductivity; θ is the volume of the soil; Lf is the latent heat of fusion of water per unit mass; and t is time. The subscripts p, w, i, and a represent the soil particles, water, ice, and air, respectively. Equation 34 does not contain the terms with respect to the heat transport through water flow and the mechanical work.

4.2 Hydraulic field

Recently published research that used molecular dynamics [41] indicated that the unfreezable threshold corresponds to pore diameter and the adsorption of the soil. This is further evidence for the analogy of the freezing and drying processes [2]. Therefore, Richard’s equation can be used to describe the fluid movement with a source term related to ice formation [3]:

$$C\frac{\partial h}{\partial t} = \nabla \cdot \left( {k\nabla \left( {h + {\mathbf{i}}} \right)} \right) - \frac{{\rho_{\text{i}} }}{{\rho_{\text{w}} }}\frac{{\partial \theta_{\text{i}} }}{\partial t}.$$
(36)

The simplified van Genuchten [35] equation for the soil–water characteristic curve (SWCC) can be extended to describe the relationships among the unfrozen water content, matric potential, and hydraulic conductivity:

$$C = \frac{\alpha m}{1 - m}\left( {\theta_{\text{s}} - \theta_{\text{r}} } \right)S_{\text{e}}^{1/m} \left( {1 - S_{\text{e}}^{1/m} } \right)^{m} ,$$
(37)
$$k = k_{s} S_{e} \left[ {1 - \left( {1 - S_{e}^{1/m} } \right)^{m} } \right]^{2} ,$$
(38)
$$S_{e} = \frac{{\theta_{w} - \theta_{r} }}{{\theta_{s} - \theta_{r} }} = \left( {1 + \left| {\alpha h} \right|^{{1/\left( {1 - m} \right)}} } \right)^{m} .$$
(39)

where h is the matric potential (negative water suction); C = dθw/dh is the specific moisture capacity; k and ks are the hydraulic conductivities for unsaturated and saturated soils, respectively; Se is the effective saturation; and θs and θr are the volumes of the saturated and residual water, respectively. Water density ρw = 1000 kg/m3 and ice density ρi = 931 kg/m3. i is the unit vector in the direction of gravity. α and m are the empirical parameters of the SWCC.

Assuming thermodynamic equilibrium at the ice–pore–water interface is maintained in an infinitesimal time interval, the rate of ice content can be derived using the chains rule with the generalized Clapeyron equation, as follows:

$$\begin{aligned} \frac{{\partial \theta_{\text{i}} }}{\partial t} & = - \frac{{\rho_{\text{w}} }}{{\rho_{\text{i}} }}\frac{{{\text{d}}\theta_{\text{w}} }}{{{\text{d}}h}}\frac{{{\text{d}}T}}{{{\text{d}}t}} = - \frac{{\rho_{\text{w}} }}{{\rho_{\text{i}} }}C\frac{{L_{\text{f}} }}{{gT_{0} }}\frac{{{\text{d}}T}}{{{\text{d}}t}},\quad T < T_{\text{f}} \\ T_{\text{f}} & = T_{0} + \frac{{gT_{0} }}{{L_{\text{f}} }}h, \\ \end{aligned}$$
(40)

where T0 = 273.15 K is the nominal freezing temperature and Tf is the freezing temperature for unsaturated soil.

4.3 The stress–strain field

The governing equation for the stress field is the Navier equation, which incorporates the equation of motion, strain–displacement correlation, and constitutive relationship. For a static problem, the equation of motion can be written in the general tensor format as:

$$\nabla \cdot \left( {{\mathbf{D}}\nabla {\mathbf{u}}} \right) + {\mathbf{F}} = 0,$$
(41)

where u is the displacement vector; D is the fourth-order tensor of material stiffness; and F is the body force vector.

The strain–displacement equation is:

$${\varvec{\upvarepsilon}} = \frac{1}{2}\left[ {\nabla {\mathbf{u}} + \left( {\nabla {\mathbf{u}}} \right)^{\text{T}} } \right].$$
(42)

The constitutive equation is:

$${\varvec{\upsigma}} = {\mathbf{D\varepsilon }}_{\text{el}} = {\mathbf{D}}\left( {{\varvec{\upvarepsilon}} - {\varvec{\upvarepsilon}}_{\text{vh}} - {\varvec{\upvarepsilon}}_{\text{p}} } \right),$$
(43)

where σ is the Cauchy stress tensor and ε, εel, εvh, and εp are the infinitesimal strain tensor, elastic strain tensor, volumetric plastic strain tensor by frost heave, and plastic strain tensor by stress, respectively.

4.4 General boundary condition

The general boundary conditions, which include one or more of the Dirichlet, Neumann, and Robin boundary conditions, are expressed as:

$${\mathbf{n}} \cdot \left( {c\nabla u + \zeta u - \gamma } \right) + qu = \delta - {\mathbf{h}}^{\text{T}} {\varvec{\upmu}},$$
(44)

where n is the outward normal unit vector of a boundary; u is the dependent variable of the individual field; c is a conductivity term; ζ is the conservative flux convection coefficient; γ is the source in the subdomain; q is the boundary absorption coefficient; δ is the boundary source; hT is a matrix representing the flexibility of the constraint type; and μ is the matrix of the Lagrange multiplier.

5 Model implementation

Equations 1 and 34–40 express the variables and the governed equations that frame heat–moisture–deformation coupling in frozen soil. A frost heave model that accounts for the lenticular structure of the fine-grained soil can be developed by combining Eqs. 2, 11, 12, 23–25, and 28–33. The equations must be solved numerically because of the high nonlinearity of the variables and the immediate updating of the local coordinate system for each element with respect to the temperature gradient. The coupled multiphysical model was solved by using COMSOL for the cross-platform finite element analysis of the coupled systems of partial differential equations.

6 Validation of the transversely isotropic frost heave model

A standard laboratory uniaxial frost heave test and field measurements of canals were performed to validate the TISO model for frozen soil. Moreover, a generalized ISO frost heave model was applied for comparison. For the ISO constitutive model, the dimensionless quantity, ξ, was set as 1/3, and the elastic constants were given by the general rule of mixtures:

$$\begin{aligned} E_{\text{ISO}} & = \theta_{\text{i}} E_{\text{i}} + \left( {1 - \theta_{\text{i}} } \right)E_{\text{m}} \\ \nu_{\text{ISO}} & = \theta_{\text{i}} \nu_{\text{i}} + \left( {1 - \theta_{\text{i}} } \right)\nu_{\text{m}} . \\ \end{aligned}$$
(45)

For ISO plasticity, yield stresses satisfied σys1 = σys2 = σys3 = σys0 and τy23 = τys31 = τys12 = τys0.

6.1 Simulation case I: standard uniaxial frost heave test and simulation

First, two of the experimental results presented by Penner [30] were repeated using numerical simulations. Both the length and the diameter of the soil column were 10 cm. The level of the external water supply was controlled to be aligned with the top of the cell. The temperature ramp rate was − 0.50 °C/day, and the temperature gradient was 0.33 °C/m in Test 1. In Test 2, the temperature ramp rate was − 0.50 °C/day with a temperature gradient of 0.12 °C/m. The initial temperature of the cold end was set at 0 °C. The sides of the samples were thermally insulated. Moreover, an axial load of 50 kPa was applied. The simulation parameters listed in Table 1 were set according to the experimental data and the relevant literature [26]. The dimensionless quantity, ξ, was 0.9 for anisotropic deformation in the TISO model [22]. Under such a small axial load, the plastic behavior of the soil was not considered.

Table 1 Soil index properties of test samples
Full size table

Comparisons of the experimental and simulated results are shown in Fig. 3. The TISO and ISO models indicated that the total frost heaving increases approximately linearly with time. Obviously, the TISO model agreed more accurately in both tests. The absolute errors for the TISO model were 0.64 mm in Test 1 and 0.57 mm in Test 2, whereas those for the ISO model were − 2.90 mm in Test 1 and − 1.29 mm in Test 2.

Fig. 3
figure 3

Comparison among the simulation results using the ISO and TISO models and the experimental results for total heave versus elapsed time. a Test 1, b Test 2. ISO isotropic, TISO transversely isotropic

Full size image

6.2 Simulation case II: monitoring and simulation of prototype canal frost heave

Li [17] observed a prototype of frost heave in a large U-shaped concrete canal. The test was conducted at the loess terrace on the east coast of the Qianhe River in the Fengjia mountain irrigation area in Baoji, China. The test canal flows from west by south to east by north, where sunny and shady slopes are clearly distinguished. The designed discharge is 58 m3/s. The canal depth is 5.1 m, with a 7.1-m-wide opening. The arc radius of the canal bottom is 3.2 m, and the canal is inclined against a wall at 10°. The canal featured a 10-cm thickness of shotcrete lining. Eleven points of measurement, N1–S1, were symmetrically deployed in the canal basin for observing and measuring the frozen depth and frost heave. Figure 4a illustrates the cross section of the canal and the distribution of the measurement points.

Fig. 4
figure 4

Test section of a canal in the Fengjia mountain irrigation area in Baoji, China. a Section dimensions (m) of the canal and the positions of the 11 measurement points, N1–S1. b Meshing and boundary conditions that control the finite element model

Full size image

Four rows of water inlets with an 8-cm-diameter and 1.5–4-m-deep opening were established on both sides of the test canal to supply water and soak the soil in the canal with water for two months before freezing occurred in the canal. The perched water table was observed 3.5–4 m beneath the canal after water was supplied. Observation started from November 14 to March 28 of the subsequent year when the daily average temperature was higher than 0 °C. Maximum frost depth and frost heave were observed on January 24.

The observed soil was silty clay containing 6.8–12.8% gravel, 55.2–62.6% fine particles, and 27.0–38.0% clay, indicating frost-susceptibility. The following properties were also experimentally determined: specific gravity = 2.71; liquid limit = 29.2–31.6%; plastic limit = 17.7–19.0%; cohesive force = 11–68 kPa; and the internal friction angle = 24.5°–31.5°. Regarding thermodynamic parameters, Cp = 3.09 × 106 J m3/K and λp = 1.85 W/(m K). In terms of hydraulic parameters, θs = 0.41, θr = 0.03, α = 0.21 m−1, m = 0.22, and ks = 3 × 10−6 m/s. Em = 10.2 MPa and υm = 0.44. The dimensionless quantity of frost heave for loess was ξ = 0.7 [39]. Other fundamental physical constants are listed in Table 1. The basic material parameters of concrete are: density ρc = 2300 kg/m3, thermal conductivity λc = 1.8 W/(m·K), volumetric heat capacity Cc = 2 × 106 J·m3/K, elastic modulus Ec = 25 GPa, and Poisson’s ratio υc = 0.33.

For the TISO plasticity of the frozen soil, σys1 = 0.2877 |T| + 0.5, σys2 = 0.2549 |T| + 0.5, τys12 = 0.18 |T| + 0.10, and τys31 = 0.425 |T| + 0.10. The unit of yield stress is MPa, and if the soil temperature is greater than 0 °C, only constant terms are used. It is noted that the slope terms of the yield shear stresses, τys12 and τys31, were determined from the adfreeze bond strength test and shear test, respectively. For the ISO model, the yield stresses are conventionally taken as σys0 = σys2, τys0 = τys31.

A two-dimensional thermal–moisture–deformation coupled finite element model of the test canal was developed. Analysis was performed for 72 days, which began on November 14 to January 24 of the subsequent year. Figure 4b shows the mesh and boundary condition control of the finite element model.

The origin of the global Cartesian coordinate system of the model was consistent with the circle in the bottom arc of the canal. The dimensions of the model are of 20 m in length and 13.38 m in width. According to Saint–Venant’s Principle and numerical analysis experience, uneven distribution of hydraulic, thermal, and stress fields in the canal only occurred within one to two times the area of canal sections. The left and right boundaries of the model were kept to be thermally insulated (∂T/∂n = 0) and impermeable (∂h/∂n = 0), and the left and right and bottom boundaries were set as a directional hinged constraint (∂u/∂n = 0). It was assumed that the underground water level was located 3.2 m from the bottom of the canal (y = − 7 m) and, before calculation, the water content of soil had reached equilibrium under gravity. Therefore, Richards’ equation for a steady-state solution was used to obtain the initial hydraulic field of the base soil. It was assumed that the soil had completed consolidation before testing. Therefore, the initial stress field of the model could be determined by using the gravity static analysis method. Boundaries of the initial soil temperature and bottom constant temperature were consistent with the temperature of the supplied water, that is, 13 °C. Ambient temperature was described in terms of a trigonometric function based on the measurement results of the three groups of daily average temperature as follows:

$$T_{\text{amb}} = 9.584 + 17.46{ \sin }\left( {\frac{2\pi }{365}\left( {t - 26.5} \right) - \frac{\pi }{2}} \right),$$
(46)

where Tamb (°C) denotes the daily average ambient temperature of the upper surface of the model and t(d) denotes the number of days in a cycle starting from January 1. Combined with the observation date, the equation yields t = 318 days when the observation commences and t = 389 days when frost depth and frost heave reach peak values.

Under the influence of solar radiation and a local winter monsoon, the canal test exhibited distinctive characteristics of sunny and shady slopes. To eliminate the effects of the two aspects on the hydraulic, thermal, and stress fields, the convective heat transfer coefficient (hc) along the upper surfaces was determined using an inversion thermal–moisture–deformation coupled finite element analysis through the observed frost depth data at t = 389 days. The obtained convective heat transfer coefficients (unit: W/m2 K) for each measurement point are as follows: hcN1 = 0.45, hcN2 = 0.30, hcN3 = 0.60, hcN4 = 1.10, hcN5 = 1.40, hc0 = 1.85, hcS5 = 1.50, hcS4 = 1.60, hcS3 = 1.28, hcS2 = 0.95, hcS1 = 1.45. The convective heat transfer coefficients of the ground surfaces in the north and south coast of the canal are hcNW = 0.55 and hcSE = 0.99, respectively. The inversion results, displayed in Fig. 5, indicate that the calculated frost depth and the observed results are basically consistent, with an absolute error ranging from − 3.38 cm to 2.63 cm.

Fig. 5
figure 5

Comparison of frost depth between simulation and test at t = 389 days

Full size image

Similar to the comparative analysis method adopted in Case I, the ISO model and TISO model of frozen soils were compared. The numerical results of the normal frost heave in the lining of the canal and the normal stress at the bottom of the lining at t = 389 days are shown in Fig. 6. The numerical results demonstrate that the normal frost heave deformation and the normal stress distribution considering TISO are both fundamentally consistent with actual measurements. The normal frost heave of the canal lining increased with frost depth, and the frost heave of the shady (north-facing) slope was greater than that of the sunny (south-facing) slope. Maximum frost heave appeared along the top of the canal because of a lesser lining constraint and its double-direction freezing of the lining surface and ground surface. The bottom of the canal exhibited relatively distinctive frost heave deformation because of the frost depth and water supply. At measurement point 0, the absolute error of the normal frost heave and normal stress were 0.64 mm and − 25 kPa. The normal frost heave and normal stress were 8.96 mm and − 104 kPa with respect to the actual measurement 8.32 mm and − 129 kPa, respectively.

Fig. 6
figure 6

Simulation results along the lining at t = 389 days when the frost heave reaches a maximum, and the variation in normal displacement of the measurement points with elapsed time. a Frost heave comparisons. A positive value indicates that the lining deformed toward the inner portion of the canal; otherwise, it occurred toward the base. b Normal stress comparisons at the bottom of the canal lining. A positive value indicates that the bottom of the canal lining suffered tension; otherwise, it suffered compression. c Normal displacement with respect to the TISO model versus elapsed time. TISO transversely isotropic

Full size image

The frost heave analysis result using the ISO model showed a nonpositive correlation between frost heave and frost depth at the measurement points. The U-shaped lining was subjected to symmetrical extrusion, which caused the bottom of the canal to deform toward the base soil. The normal frost heave at measurement point 0 at the bottom of the canal was − 3.52 mm, and the actual measurement was 8.32 mm, indicating an absolute error of − 11.84 mm.

As shown in Fig. 6c, the amount of frost heave was positive and was primarily related to the depth of freezing. However, due to asymmetric frost heaving of the shady slope and sunny slope, the canal lining tended to tilt toward the sunny slope. At the beginning of the frost heave (t = 357–373 days), N4 showed a distinct normal displacement toward the foundation because of the hysteresis of freezing at the sunny slope. After that, the lining played a certain role of constraint in the frost heaving process of the foundation soil, especially at the canal bottom. Also, the normal displacements of the measured points at shady slope (N5–S2) exhibited a stepped increase with freezing time.

The effective plastic strain distribution calculated with the elastoplastic model at t = 389 days is shown in Fig. 7. It can be seen that the ground frost within a 70 cm depth received various degrees of yielding due to the lining confining. The effective plastic strain occurred in the range of N3–S2. The maximum of the effective plastic strain was approximately 0.019 and occurred at approximately point S5.

Fig. 7
figure 7

Effective plastic strain distribution obtained by using the elastoplastic constitutive in the TISO model within 70 cm of soil depth at t = 389 days. TISO transversely isotropic

Full size image

7 Discussion

7.1 Equivalence of micro-ice lenses and soil

In the derivation of TISO elastic constants, the freezing of microstructures into ice lenses was considered. However, it was found that an accurate approach needs to be developed to distinguish distribution problems concerning pore ice and ice lens volume. The mixing rule model should be adopted for individual elastic constants estimations of the soil layers with dispersive pore ice. Hence, series models and parallel models should be further employed for finding the equivalent elastic behavior of frozen soil that contains ice lenses and soil layers. However, in freezing fine-grained soils, the soil layers between ice lenses showed dewatering rather than a constant saturation during decreasing temperatures. Therefore, this influence was not considered in the present study.

7.2 Quantitative analysis of transversely isotropic elastic constants and ice contents

A comparison of E1 in Eq. 23 and the mixing rule calculation, Eq. 45, shows that the first two terms in both equations are identical, but Eq. 23 contains the correction term of the elastic constants and volumes for the ice and soil. A dimensional analysis revealed that this value was slightly smaller than the square (1%) of the difference of the soil–ice Poisson’s ratios. The removal of this correction term led to a negligible error when compared with the elastic modulus (MPa) of the soil matrix, Em. Similarly, a comparison of υ12 in Eqs. 24 and 45 showed that the equivalent Poisson’s ratio obtained using the mixing rule equation was identical to that of the first two terms in the υ12 equation. The difference was also the square of υ. However, compared with υ, it bore the same weight and thus should not be omitted. In contrast, a greater difference of the elastic constants was obtained in other main directions.

Figure 8 illustrates the relationship among the change of ice content, θi, the TISO elastic modulus, and Poisson’s ratio when Em = 46 MPa, υm = 0.4, and Ei = 6000 MPa, υm = 0.3. Figure 8a shows that the curve of the elastic modulus, E1, is parallel to the direction of the ice lenses, which increases linearly with the ice content. As the correction term of E1 is omissible, based on the mixing rule model, there is an overlap between E1 and the equivalent elastic modulus EISO. The curve of E2, which is vertical to the direction of ice lenses, performs as an S-shaped function and increases with the ice content, θi. In this case, the two inflection points appeared at θi = 0.05 and θi = 0.8, respectively. Within the range of θi from 0.05 to 0.8, E1 ≫ E2 and had no significant change.

Fig. 8
figure 8

Relations between the elastic constants of frozen soils and ice content described using the ISO and TISO models with a soil elastic modulus Em = 46 MPa, Poisson’s ratio υm = 0.4, and ice with Ei = 6000 MPa, υi = 0.4. a A comparison of the elastic moduli with increasing ice content. b A comparison of Poisson’s Ratio with increasing ice content. ISO isotropic, TISO transversely isotropic

Full size image

Similarly, Fig. 8b shows that Poisson’s ratio, υ12, first increases and then slowly decreases to the value of the ice Poisson’s ratio, υi. υ21 quickly dropped to υi when θi < 0.2 and thereafter showed no significant change as θi increases. The variation magnitude for υ13 is maximized for increasing θi, and υ13 quickly decreases when θi < 0.2. When θi is between 0.2 and 0.8, υ13 exhibits minimal change (almost a straight line), quickly regressing to υi when θi > 0.8. The equivalent Poisson’s ratio, υISO, based on the mixing rule model decreases linearly to the Poisson’s ratio of ice as θi increases. When θi = 0.1, υISO has maximum difference with the Poisson’s ratio obtained in the present study. The maximum difference between υISO and υ13 is approximately 0.32.

When θi is between 0.2 and 0.8, Poisson’s ratio of the ice lens, υ13, reflects a horizontal straight line. This result matches well with the results given by Arenson et al. [1] who reported that frozen soil can be considered a TISO material when ice makes up 10–90% of the soil. Therefore, additional verification is required to judge the feasibility of determining the value of υ13 using Eq. 25 to make sure the frost-susceptible soil easily forms ice lenses, or to analyze the problem of distinguishing pore ice and ice lenses by using the initial and final points of the horizontal section of a curve.

7.3 Interaction between ice and soil layers

The assumption of an ideal bond between the ice layers and soil layers is essential for the derivation of elastic constants in this study. In fact, the clay characteristics of frozen soil (e.g., unfrozen water content) and temperature directly influence the bond and coordinated deformation. Hence, when frozen soil temperature and unfrozen water content are high during the initial freezing stage, the ice lens is affected by unfrozen water film, which may result in a minimal nonsynchronous slip or cause soil particles to move. This slip process deviates somewhat from this assumption. However, when frozen soil is constrained by overburden pressure and ice crystals of the freezing fringe, the effect of this slight deviation is negligible.

When the stress of the frozen soil is greater than the yield stress, failure (bend, fold, and fracture) of ice lenses in the soil easily occurs, with a possible slip between ice lenses and the soil bodies. In fact, the mechanism of the interaction between ice lenses and soil is complex. Therefore, a reasonable approach should be developed to examine and establish a standard thermal–moisture–deformation coupled transversely ISO nonlinear constitutive model.

7.4 Frost heave partition coefficient ξ

Uneven frost heaving of frozen soil involves complex coupling of multiple physical fields and is closely related to the particle grade of soil and the interaction between soil particles and ice crystals. The mechanism of uneven frost heaving has never been reported. Estimations of frost heave partition coefficients can only be obtained by using specific tests or parameter inversion.

The calculation results of case I show that the frost heave partition coefficient is not sensitive to the cooling rate for a specified type of soil. In case II, the frost heave partition coefficient can be approximated as a constant in the two-dimensional frost heaving analysis of complex temperature boundaries.

8 Conclusions

Based on the natural frost heaving deformation process and mesoscopic composition of fine–grained soil, with the assumption of an ideal bond between ice and soil, a typical micro-soil–ice series model and parallel model for frozen soil with ice lenses were established in this study by using deformation coordination, force equilibrium, and Hooke’s law. Subsequently, the TISO constitutive model of frozen soil was presented, and the expressions for the five independent elastic constants were derived on the basis of mechanics of composites. The following conclusions were drawn.

First, an improved theoretical method for solving TISO elastic constants was improved in this study. It overcame the disadvantage in the mixing rule theory, which can only approximately describe the elastic constant in symmetry in the plane (planes 1–3), but fails to reflect the other four elastic constants of transversely ISO frozen soil.

Second, a preliminary elastoplastic model for TISO frozen soil was established based on Hill plastic theory.

Third, a TISO frost heave model and numerical method for frozen soils were established based on the thermal–moisture–deformation coupling theory. Synchronous coupling and uniformity of the TISO constitutive model and the TISO frost heave model were realized for frozen soils with alternately laminated ice lenses structures. This approach overcame the drawbacks of previously used models, such as the rigid ice and segregated ice models, which are only applicable for one-dimensional freezing of saturated soil. Therefore, this approach bridged the research gap of the hydraulic frost heave model, which cannot consider the TISO mesoscopic composition of frozen soil and the uneven frost heave evolution.

Fourth, a standard indoor ramped freezing test together with the freezing process of a prototype canal freezing process was simulated using the thermal–moisture–deformation coupling finite element method. Results indicated that the TISO constitutive model for freezing fine-grained soil was more accurate than the ISO constitutive model. The analytical expressions of the elastic constants for frozen soil with ice lenses derived in this study are essential for deeply understanding the interaction between freezing fine-grained soils and buildings in a complex environment and with complex boundaries.

Finally, the stress equations for ice lenses and soil pressure were proposed for estimating ice pressure in frozen soil.