Abstract
The deformation property of marine clay under a heat source has received considerable attention in the geotechnical literature. In this paper, a three-parameter fractional order derivative model is introduced into the thermo-hydro-mechanical coupling governing equations with thermal filtration and thermo-osmosis to simulate viscoelastic characteristics of marine clay. The excess pore pressure, temperature increment, and displacement of marine clay are derived by using the Laplace transform method, and the semianalytical solution for the one-dimensional thermal consolidation in the time domain is derived by using a numerical inversion of the inverse Laplace transform. The influence of the order of the fractional derivative, material parameters, and phenomenological coefficient on thermal consolidation is investigated based on the present solutions. It is shown that the influence of the fractional derivative parameter on the excess pore pressure and displacement of marine clay depends on the properties of soil mass, and the temperature increment has an obvious effect on the thermal filtration and thermo-osmosis process.
1. Introduction
Over the past decades, more and more engineering activities have caused the change in the temperature field of the surrounding soil layers, which inevitably affects the physical and mechanical properties of the strata. These engineering activities involve a broad range of civil engineering topics such as deep geological disposal of radioactive waste [1], deep drilling and excavation [2, 3], extraction of geothermal energy [4–6], energy piles [7, 8], ground improvement using prefabricated vertical thermal drain [9–11], oil and gas pipelines [12], and frictional heating-induced large-scale landslides [13]. This huge engineering demand has stimulated scholars to pay their attention on the thermo-hydro-mechanical coupling theory of porous media, especially the deformation properties of marine clay under a heat source [14].
Biot originally proposed the constitutive equations considering the thermal effect and established the theoretical framework for the thermodynamics of saturated porous media [15]. Following Biot’s research, a number of scholars have investigated the thermal consolidation of elastic saturated porous media by considering the influence of the thermal effect. Booker and Savvidou studied the thermal consolidation of saturated soil under the action of the deep buried spherical and point heat source [16, 17]. Bai [18] derived an analytical solution for the wave response of porous materials under cyclic thermal load [18]. Smith and Booker [19] established the linear thermoelastic theory of homogeneous isotropic materials by using the Green function method [19]. Lu et al. [20] investigated the thermal-mechanical coupling response of saturated porous media under simple harmonic heat and load based on the generalized thermoelastic theory [20]. Ai and Wang [21] studied the axisymmetric thermal consolidation of layered elastic saturated porous media under the action of the heat source [21]. Furthermore, Ai et al. [22] studied the axisymmetric thermal consolidation of layered transversely isotropic porous media [22]. Wen et al. [23] studied the thermal dynamic response of a lined circular tunnel in saturated elastic porous media with the help of a fractional order derivative [23].
Although the above-mentioned scholars have greatly promoted the development of the thermo-hydro-mechanical coupling theory of porous media, they neglect thermal filtration (the influence of the excess pore pressure gradient on heat flux) and thermo-osmosis (the influence of the temperature gradient on water flux). In reality, by conducting an experiment on thermo-osmosis of water through kaolinite, Srivastava and Avasthi [24] found that the water flux associated with thermo-osmosis in compacted kaolinite can reach 10−8 m/s under a temperature gradient of 20°C/m [24]. Considering the influence of thermal filtration and thermo-osmosis, Zhou et al. [25] presented a fully coupled thermo-hydro-mechanical model [25]. Considering the influence of seepage velocity on thermal diffusion and the temperature gradient on seepage velocity, Liu et al. [26] established a fully dynamic coupled thermo-hydro-mechanical model and investigated the influence of the thermal and water permeability coefficient of saturated porous media around cylindrical holes and spherical cavities on their thermodynamic response, respectively [26, 27]. Yang et al. [25] proposed a refined mathematical model to study the coupled effect of thermo-osmosis in saturated porous media [28]. All these references utilized the elastic theory to simulate stress-strain relationships in saturated porous media, which may not be applicable to the marine clay covering two-thirds of the Earth.
With the continuous development of marine development all over the world, more and more attention has been paid to the deformation characteristics of marine clay. Liu et al. [29] studied the one-dimensional consolidation of viscoelastic marine soft soil under load varying with depth and time [29]. Liu et al. [14] investigated the effect of viscosity on the one-dimensional thermal consolidation of marine soft soil [14]. Wang and Wang [30] studied the rheology and thermal consolidation of layered saturated soft soil under force and thermal load [30]. Although these works have promoted the understanding of the deformation characteristics of marine clay, the existing research is still far from being applied to the engineering applications in marine clay due to its complex rheological properties. In order to reasonably consider the rheological characteristics of soil, many scholars introduced the fractional constitutive model to study the consolidation characteristics of soil [31–35]. However, to the authors’ knowledge, there is no report on the thermal consolidation of viscoelastic marine soil with a fractional order derivative.
In light of the above, the objective of this paper is to investigate the one-dimensional thermal consolidation of viscoelastic marine clay with the fractional order derivative. By introducing a three-parameter fractional order derivative model to consider thermal filtration and thermo-osmosis, the solutions of excess pore pressure, temperature increment, and displacement of viscoelastic marine clay are obtained by using the Laplace transform method and its numerical inverse transform. Based on the present solutions, the influence of the order of the fractional derivative, material parameters, and phenomenological coefficient on the thermal consolidation of viscoelastic marine clay is investigated.
2. Mathematical Modeling
Considering the effects of thermal filtration and thermo-osmosis, a fully coupled thermo-hydro-mechanical model was proposed by Zhou et al. [25]. On this basis, a three-parameter fractional order derivative model, consisting of two springs and one dashpot in parallel, is introduced into the thermo-hydro-mechanical coupling governing equations. The constitutive relationship of the fractional order derivative model is expressed as [36]Here, denotes the vertical effective stress; and are Lame constants; represents the depth below the ground surface; t is time; , and are material parameters; and is the fractional order parameter which can be obtained by the triaxial test and parameter inversion method, and ; is the displacement in the z direction. denotes order Riemann–Liouville fractional derivative, and it is defined as [36]Here, is the Gamma function.
As shown in Figure 1, both the top and bottom surfaces are set as permeable. Instantaneous compressive stress q0 and temperature T1 are applied at the top surface, and a constant temperature T0 is maintained at the bottom surface.

The viscoelastic characteristics of marine clay are simulated by the three-parameter fractional order derivative model, the equilibrium equation is given by [25]Here, depends on the compressibility of the soil grains; represents the drained bulk modulus of the soil medium; denotes the bulk modulus of the soil grains; is the coefficient of volumetric expansion of the soil medium and T is the temperature increment; and is excess pore pressure.
Substituting equation (1) into equation (3) yields
The fluid mass balance is governed byHere, k denotes the coefficient of permeability; is a phenomenological coefficient associated with the influence of the thermal gradient on the water flux (thermo-osmosis); the parameters c1, c2, and c3 are expressed as , , and , respectively; n is the porosity; and are volumetric thermal expansion coefficients of the pore water and soil grains, respectively; and is the bulk modulus of pore water.
The thermal energy balance equation is given aswhere is the volumetric specific heat of the soil medium; ; is the density of soil grains; is the density of water; and denote the gravimetric specific heats of soil grains and the pore water, respectively; and is the thermal conductivity of the soil medium, in which and are the thermal conductivities of the soil grain and the pore water.
The initial conditions of the problem are as follows:
The boundary conditions of the problem are as follows:
3. Solutions to Governing Equations
Laplace transform is introduced to solve governing equations (4)–(6), and the following Laplace transform equation is defined:
Applying Laplace transforms to equations (4)–(6), one can obtain transformed governing equations asHere, s is the Laplace transform parameter; , , and .
The Laplace transform of equations (8)–(12) with respect to the time variable t is derived as
Substituting equation (15) into equation (16) yieldsHere, ; ; ; and .
Substituting equation (15) into equation (14) giveswhere ; ; ; ; and .
Further transformation of equations (22) and (23) leads towhere ; ; and .
Substituting equation (24) into equation (22) yieldswhere ; ; and .
The general solution of equation (25) is derived aswhere A1, A2, B1, B2 are undetermined coefficients.;
Substituting equation (26) into equation (24) gives
Substituting equation (26) into equations (17), (18), (20), and (21), the undermined coefficients can be obtained aswhere , , , and .
Then, equation (26) can be rewritten aswhere, and .
Substituting equations (26) and (27) into equation (14) yieldswhere .
Integrating both sides of equation (30) from zero to infinity gives
Substituting equation (31) into equation (19), the coefficient D1 can be obtained as
Through the above derivation, the analytical expressions of excess pore pressure, temperature increment, and displacement can be obtained accordingly. However, it is difficult to directly obtain the analytical solutions in the transformed space for some formulas are difficult to integrate. Therefore, it is necessary to utilize the numerical method to program and analyze the solution process in this paper. Currently, there exist many Laplace inverse transform methods, among which the Crump’s numerical inversion results of Laplace transform is the most accurate [10].
Assuming that F (s) is the Laplace transform of function F (t), Crump’s inversion algorithm of the Laplace inverse transform can be written aswhere . If , the error .
4. Numerical Results and Discussion
Having illustrated and reviewed the obtained theoretical results in the preceding section, we now present some numerical results. Following the works of Liu et al. [14] and Zhou et al. [25], the parameters of marine clay in the calculation can be set as Pa, , , m2/s/°C, Pa, m5/J/s, J/m/s/°C, J/m/s/°C, Pa, Pa, kg/m3, kg/m3, J/kg/°C, J/kg/°C, °C−1, °C−1, , Pa, °C, and °C.
4.1. Comparative Analysis
The Kelvin viscoelastic model is combined with the thermo-hydro-mechanical coupling governing equations for marine clay, and the influence of viscosity on one-dimensional thermal consolidation was investigated by Liu et al. [14]. In addition, Zhou et al. [25] proposed a fully coupled thermo-hydro-mechanical model with thermal filtration and thermo-osmosis. In this paper, for the purpose of verifying and presenting the differences of the three cases of the fractional derivative viscoelastic model (FDVM), Kelvin viscoelastic model (KVM), and elastic model (EM), the order of the fractional derivative is assumed to be zero for the comparison of the FDVM with the EM, and the material parameter is assumed to be zero, and is equal to 1 for the comparison of the FDVM with the KVM. Figures 2–4 illustrate the distribution of excess pore pressure, displacement, and temperature increment against depth among solutions of Zhou et al. [25], Liu et al. [14], and present work. It can be found that the KVM produces the largest excess pore pressure, followed by the results of the FDVM and EM. At depth above 0.6 m, the displacement value produced by the EM is the largest, followed by those of the FDVM and KVM. At larger depth, remarkable increase of the displacement value could be observed for all the three models. For the distribution of temperature, the three models produce similar results.



4.2. Influence of the Order of the Fractional Derivative
Figures 5 and 6 show the influence of the parameter on the development of excess pore pressure when the values of are 3.0 and 0.5, respectively. When , the excess pore pressure increases with the increase of before the peak value. After the peak, the increase of seems to accelerate the dissipation of the excess pore pressure. When , an opposite trend could be observed. With increasing , the peak value of excess pore pressure shows a decreasing trend, and the after-peak dissipation rate is smaller. It can be seen that the influence of the fractional derivative parameter on the development of excess pore pressure depends on the characteristics of the soil mass.


Figures 7 and 8 show the influence of the parameter on the displacement development under the same conditions. Similar trends to that of the excess pore pressure can be observed. Another finding is that the influence of is more obvious at . It can be seen that the soil property would influence the impact of the fractional derivative parameter on the development of displacement.


4.3. Influence of Material Parameters
Figures 9 and 10 depict the influence of variance on the development of excess pore pressure and displacement with . With increasing , both the peak value of the excess pore value and the after-peak dissipation rate are much larger, while the net maximum displacement value decreases. The influence of on the trend of displacement development is less obvious.


The influence of embedment depth on the development of excess pore pressure and displacement is presented in Figures 11–14. As shown in Figures 11 and 12, the magnitude of excess pore pressure at is much larger than that at . On the other hand, the development trend of excess pore pressure at different depths seems to be uninfluenced with the variance of , but the peak excess pore pressure decreases with depth. As depicted in Figures 13 and 14, similar trends could be observed for the development of the displacement value. With increasing depth, the displacement value is much smaller.




4.4. Influence of the Phenomenological Coefficient
The influence of the phenomenological coefficient on the development of excess pore pressure, displacement, and temperature increment is presented in Figures 15–17. It is clear that both the excess pore pressure and displacement increase with increasing . On the other hand, the variance of seems to have no impact on the development of temperature increment. This parametric analysis indicates that the thermal gradient has an obvious effect on the thermal filtration and thermo-osmosis process.



5. Conclusion
In this paper, a three-parameter fractional order derivative model is proposed based on the improved thermo-hydro-mechanical coupling theory by accounting for the viscoelastic behavior of marine clay. One-dimensional thermal consolidation of viscoelastic marine clay is analyzed, considering the thermal filtration and thermo-osmosis process. Solutions of excess pore pressure, temperature increment, and displacement are obtained by using the Laplace transform method and its numerical inverse transform. The influence of the order of the fractional derivative, material parameters, and phenomenological coefficient on the characteristics of thermal consolidation of marine clay is investigated. The performance of the proposed fractional derivative viscoelastic model (FDVM), Kelvin viscoelastic model (KVM), and elastic model (EM) is compared. The calculated excess pore pressure values from the three models are quite different, while the development of temperature increment seems to be uninfluenced by the selected model. For the development of the displacement value, there is an obvious turn of the displacement value for different models at the depth of 0.6 m. Further parametric analysis indicates that the influence of the fractional derivative parameter on the development of excess pore pressure and displacement depends on the properties of the soil mass, and the temperature increment has an obvious effect on the thermal filtration and thermo-osmosis process.
Data Availability
Data are available on request from the authors.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (Grant nos. 52108347, 52178371), the Outstanding Youth Project of the Natural Science Foundation of Zhejiang Province (Grant no. LR21E080005), and the Exploring Youth Project of Zhejiang Natural Science Foundation (Grant no. LQ22E080010). The China Postdoctoral Science Foundation Funded Project (Grant no. 2020M673093) and the Construction Research Founds of Department of Housing and Urban-Rural Development of Zhejiang Province (Grant no. 2021K256), and the funds from the Engineering Research Center of Rock-Soil Drilling & Excavation and Protection, Ministry of Education (Grant No.202203) are also acknowledged.