Dynamic response analysis of layered saturated frozen soil foundation subjected to moving loads

Abstract

Based on the theory of solid porous media with pores, the dynamic response of layered saturated frozen soil under uniform moving load is studied. Firstly, the governing equation of saturated frozen soil is established. Then, the Fourier integral transform is used to decouple the governing equation of saturated frozen soil, and the general solution of the potential function is obtained, and the corresponding force and displacement of each layer of saturated frozen soil are derived. Finally, using the transfer matrix method, combined with the continuous condition, the semi-analytical solution of the layered saturated frozen soil medium in the frequency domain under the surface permeable condition is derived. After verifying the accuracy of the solution by comparing the model degradation with the existing literature, the effects of soil shear modulus, load moving speed, temperature, contact parameter and frequency on the dynamic response are analyzed in detail. The results show that the order of soft and hard soil layers is of great significance to accurately evaluate the dynamic response of saturated frozen soil.

Keywords: Saturated frozen soil; Fourier integral transformation; Transfer matrix method; Dynamic response

Graphical Abstract

1 INTRODUCTION

Frozen soil is prevalent across various regions worldwide, predominantly concentrated in Europe and northern Asia. With the rapid development of economy, many infrastructure projects such as bridges, highways and railways have been built, under construction and proposed in permafrost regions. With the continuous improvement of the basic transportation network, the vibration problem of frozen soil foundations caused by high-speed railway and automobile traffic is becoming more and more serious. Frozen soil changes the constitutive relationship and general mechanical properties of soil due to the presence of pore ice phase, making it different from unsaturated soil and saturated soil. Therefore, in order to ensure the stability and safety of the foundation in frozen soil area, it is urgent to study the dynamic response of frozen soil under moving load.

In an early study of the dynamic response of foundations under moving loads, the site soil medium is often simplified as either an elastic or viscoelastic medium (Eason, 1965; Yang et al., 2003; Sheng et al., 2004; Beskou and Theodorakopoulos, 2011;). Due to the sedimentation of geology, the soil often presents layered characteristics. Following Thomson's pioneering (Thomson, 1950) use of the transfer matrix method (TMM) to address wave propagation in layered elastic media, many scholars have carried out theoretical research on the dynamic response of layered foundation (Li et al., 2015; Dong et al., 2021; Fan et al., 2022). Considering the influence of liquid in soil on dynamic response, Biot (1956a, 1956b, 1962) first simulated geotechnical materials as saturated media for research and proposed the wave equation of saturated media. Subsequently, more and more scholars have investigated the dynamic response of saturated ground under moving loads based on Biot's theory (Theodorakopoulos et al., 2004; Zhang et al., 2014b). Ai and Ren (2016) used the moving coordinate system to transform the dynamic problem of the moving load into a static problem, and used the analytical layer element method to study the plane strain response of the transversely isotropic saturated soil foundation surface with vertical moving load. Ba et al. (2017) used the stiffness matrix method to obtain the analytical solution of the harmonic load acting on the surface or interior of the two-dimensional foundation. Liu et al. (2022) studied the dynamic response of two-dimensional saturated soil foundation with moving strip load by using Fourier series and dual variable and position method. In recent years, in order to make the model more practical, more and more scholars have extended the foundation dynamic response problem to three dimensions. Li et al. (2020) developed a more practical transversely isotropic foundation model, and obtained the analytical solution of the foundation under moving load or pore water stress by using Fourier integral transform and stiffness matrix method. Ai et al. (2021) and Yang and Jia (2023) simplified the vehicle load as a vertical rectangular moving load, studied the dynamic response of the porous elastic medium surface with a vertical rectangular moving load, and obtained the exact solution of the dynamic response of the porous elastic medium surface with a vertical rectangular moving load.

However, the geotechnical medium in nature is partially saturated, that is, unsaturated soil, which consists of three phases and the influence of the presence of the pore gas phase on the dynamic response of the soil body cannot be disregarded (Steeb et al., 2014; Zhang et al., 2014a; Wang et al., 2015). Ye and Ai (2021) used the double Fourier transform to study the dynamic response of layered unsaturated soil under moving harmonic loads and compared it with the results of saturated porous elastic media. Tang et al., (2020) investigated the dynamic response of the coupled system of pavement and unsaturated soil foundation under traffic load. Ma et al. (2023a) explored the influence of surface soil softness and stiffness on the dynamic response of layered unsaturated foundations under moving loads using the transfer and reflection matrix method and found that the arrangement order of the soil layers had a notable impact on the surface displacement.

By combing the above-mentioned foundation dynamic response work, we can find that the above-mentioned foundation dynamic response problems often simplify the site soil medium into single-phase elastic soil, saturated soil and unsaturated soil. However, in regions with high latitudes and high altitudes, the predominant soil type is frozen soil. The dynamic response of frozen soil differs significantly from that of unfrozen soil due to the presence of the ice phase. Leclaire (1994), Leclaire et al. (1998) ignored the interaction between pore ice and porous skeleton in saturated frozen soil, established the theory of wave propagation in frozen soil, namely the LCA model, and verified the rationality of the theory by experiments. Carcione and Seriani (2001), Carcione et al. (2000, 2003) modified the saturated frozen soil model proposed by Leclaire by further considering the coupling effect between soil particle skeleton and pore ice. At present, the research on saturated frozen soil is mostly focused on wave propagation (Jiang et al., 2023a, 2023b) and the dynamic response of piles (Li et al., 2019; Cao et al., 2019; Chen et al., 2023). The layered characteristics of foundation soil are not considered. Few literatures study the dynamic response of layered saturated frozen soil foundation under moving harmonic load based on the governing equation of saturated frozen soil. Therefore, it is undoubtedly of great theoretical and practical value to study the dynamic response of layered saturated frozen soil foundation for engineering construction in frozen soil area.

In this paper, the frozen soil medium is simplified as saturated frozen soil, and the transfer matrix method is used to solve the dynamic response problem of layered saturated frozen soil under moving harmonic load. Firstly, the dynamic governing equation of 2D saturated frozen soil foundation under moving load is established. By introducing the potential function and using the Helmholtz principle, and then through the coordinate transformation and Fourier transformation, the control equation is decoupled, and the general solution of the displacement and stress of each layer expressed by the potential function is obtained. Then the analytical solution of stratified saturated frozen soil foundation subjected to moving loads in the frequency domain is derived by using the transfer matrix method and combining it with the continuity condition. Finally, the dynamic response of each phase in the time domain is obtained by Fast Fourier Transform (FFT). After verifying the model degradation and comparing it with existing literature solutions, a comprehensive analysis is conducted on the effects of soil shear modulus, load moving speed, temperature, contact parameters, and load frequency on the dynamic response.

2 GOVERNING EQUATIONS AND BOUNDARY CONDITIONS

The analysis model of layered saturated frozen soil foundation under moving harmonic load is shown in Figure 1. The frozen soil foundation is simplified to saturated frozen soil foundation for the study. Adopting the definitional approach of Zhou and Lai (2011), ice is considered to form in pores and coexist with liquid water in the pores, and assumes that each layer of the saturated frozen soil foundation as homogeneous isotropic frozen saturated porous medium consisting of the soil particulate phase, the pore ice phase, and the pore water phase. A Cartesian coordinate system is established with the horizontal direction as the $x_1(x)$ axis and the vertical direction as the $x_3(z)$ axis, where is the fixed coordinate system and the origin of the coordinates are placed on the surface of the foundation, i.e., z=0 on the surface of the foundation, and $(x,z)$ is the moving coordinate system. A layer of saturated frozen soil with a thickness H covers the impermeable bedrock. The foundation surface is under the action of a moving-line homogeneous simple harmonic load F(x,t), assuming that the moving load has a frequency of ω, a length of 2l, and an amplitude of q, $q=q_0/2l$ , and is moving uniformly with a velocity $c$ along the positive direction of the x-axis. For N-layer saturated frozen soil foundations, water is assumed to be permeable at the surface and impermeable at the underlying bedrock, and the interfaces between the layers are horizontal and well-bonded. In this paper, the parameters of the solid, liquid, and ice phases in saturated frozen soil are denoted by the upper or lower subscripts "$S$", "$F$" and "$I$", respectively.

Figure 1
A model of layered saturated frozen soil with impermeable rigid bedrock subjected to moving loads

2.1 Basic equations

The constitutive equations of the saturated frozen soil medium in the plane Cartesian coordinate system are (Carcione et al., 2000, 2003; Carcione and Seriani, 2001):

$$\begin{array}{l}{\sigma }_{ij}^S=(K_1{\theta }_S+C_{12}{\theta }_F+C_{13}{\theta }_I){\delta }_{ij}+2{\mu }_{11}{d}_{ij}^S+{\mu }_{13}{d}_{ij}^I\\ {\sigma }^F=C_{12}{\theta }_S+K_2{\theta }_F+C_{23}{\theta }_I\\ {\sigma }_{ij}^I=(C_{13}{\theta }_S+C_{23}{\theta }_F+K_3{\theta }_I){\delta }_{ij}+2{\mu }_{33}{d}_{ij}^I+{\mu }_{13}{d}_{ij}^S\end{array}$$ (1)

in which:

$${\theta }_a=u_{i,i}^a,d_{ij}^a={\epsilon }_{ij}^a-{\theta }_a{\delta }_{ij}/3,{\epsilon }_{ij}^a=(u_{i,j}^a+u_{j,i}^a)/2$$ (2)

where ${\sigma }_{ij}^S$, ${\sigma }^F$, ${\sigma }_{ij}^I$ denote the soil corresponding force, the stress borne by pore water, and ice corresponding force, respectively; $C_{12}$, $C_{13}$and $C_{23}$ are the bulk modulus between the three phases; $K_1$, $K_2$, $K_3$ are the bulk modulus of elasticity; ${\mu }_{11}$, ${\mu }_{13}$, ${\mu }_{33}$are the stiffness parameters; ${\delta }_{ij}$ is the Kroenke symbol; ${\theta }_a$, $d_{ij}^a$, ${\epsilon }_{ij}^a$ denote the a (a=S, F, I) phase body strain, partial strain, and strain, respectively; $u^a$denotes the displacement of the a (a=S, F, I) phase. ${\mu }_{11}$, ${\mu }_{13}$, ${\mu }_{33}$, $C_{12}$, $C_{13}$, $C_{23}$, $K_1$, $K_2$, and $K_3$ are shown in ^{Appendix 1} for details.

The motion equation based on solid porous media with pores can be expressed as follows (Carcione et al., 2000, 2003, Carcione and Seriani, 2001):

$$\begin{array}{l}{\sigma }_{ij,j}^S={\rho }_{11}{\ddot{u}}_i^S+{\rho }_{12}{\ddot{u}}_i^F+{\rho }_{13}{\ddot{u}}_i^I+b_{12}\text{(}{\dot{u}}_i^S-{\dot{u}}_i^F\text{)}+b_{13}\text{(}{\dot{u}}_i^S-{\dot{u}}_i^I\text{)}\\ {\sigma }_{,i}^F={\rho }_{12}{\ddot{u}}_i^S+{\rho }_{22}{\ddot{u}}_i^F+{\rho }_{23}{\ddot{u}}_i^I+b_{12}\text{(}{\dot{u}}_i^F-{\dot{u}}_i^S\text{)}+b_{23}\text{(}{\dot{u}}_i^F-{\dot{u}}_i^I\text{)}\\ {\sigma }_{ij,j}^I={\rho }_{13}{\ddot{u}}_i^S+{\rho }_{23}{\ddot{u}}_i^F+{\rho }_{33}{\ddot{u}}_i^I+b_{13}\text{(}{\dot{u}}_i^I-{\dot{u}}_i^S\text{)}+b_{23}\text{(}{\dot{u}}_i^I-{\dot{u}}_i^F\text{)}\end{array}$$ (3)

where ${\rho }_{ij}\left(i,j=\mathrm{1,}\mathrm{2,}3\right)$ are the coupling inertia coefficient between each phase;$b_{12}$, $b_{23}$, $b_{13}$ are viscosity parameters; The specific values of ${\rho }_{ij}$, $b_{12}$, $b_{23}$ and $b_{13}$ are detailed in ^{Appendix 1}.

Combining Eq. (1), Eq. (2), and Eq. (3), we can derive the governing equations for saturated frozen soil as follows:

$$\begin{array}{l}{\rho }_{11}{\ddot{u}}^S+{\rho }_{12}{\ddot{u}}^F+{\rho }_{13}{\ddot{u}}^I=R_{11}\nabla \left(\nabla \cdot u^S\right)+R_{12}\nabla \left(\nabla \cdot u^F\right)+R_{13}\nabla \left(\nabla \cdot u^I\right)-\\ {\mu }_{11}\nabla \times \nabla \times u^S-{\mu }_{13}\nabla \times \nabla \times u^I-\left(b_{12}+b_{13}\right){\dot{u}}^S+b_{12}{\dot{u}}^F+b_{13}{\dot{u}}^I\end{array}$$ (4a)

$$\begin{array}{l}{\rho }_{12}{\ddot{u}}^S+{\rho }_{22}{\ddot{u}}^F+{\rho }_{23}{\ddot{u}}^I=R_{12}\nabla \left(\nabla \cdot u^S\right)+R_{22}\nabla \left(\nabla \cdot u^F\right)+\\ R_{23}\nabla \left(\nabla \cdot u^I\right)+b_{12}{\dot{u}}^S-\left(b_{12}+b_{23}\right){\dot{u}}^F+b_{23}{\dot{u}}^I\end{array}$$ (4b)

$$\begin{array}{l}{\rho }_{13}{\ddot{u}}^S+{\rho }_{23}{\ddot{u}}^F+{\rho }_{33}{\ddot{u}}^I=R_{13}\nabla \left(\nabla \cdot u^S\right)+R_{23}\nabla \left(\nabla \cdot u^F\right)+R_{33}\nabla \left(\nabla \cdot u^I\right)\\ -{\mu }_{13}\nabla \times \nabla \times u^S-{\mu }_{33}\nabla \times \nabla \times u^I+b_{13}{\dot{u}}^S+b_{23}{\dot{u}}^F-\left(b_{13}+b_{23}\right){\dot{u}}^I\end{array}$$ (4c)

where $R_{11}=K_1+\frac{4}{3}{\mu }_{11}$, $R_{22}=K_2$, $R_{33}=K_3+\frac{4}{3}{\mu }_{33}$, $R_{12}=C_{12}$, $R_{23}=C_{23}$, $R_{13}=C_{13}-\frac{2}{3}{\mu }_{13}$(Qiu et al, 2018, Qiu, 2019); $\nabla \cdot$ is scalar product (point product), $\nabla \times$ is vector product (fork product); $R_{ij}\left(i,j=\mathrm{1,}\mathrm{2,}3\right)$ are detailed in ^{Appendix 1}.

2.2 Boundary and continuity conditions

Assuming good bonding between any adjacent layers in saturated frozen soil, for the saturated frozen soil layered system, consideration of the coordination conditions between the layers can be obtained:

$$\begin{array}{l}{\sigma }_{zz}^S(x,z_{k+},t)={\sigma }_{zz}^S(x,z_{k-},t);{\sigma }_{xz}^S(x,z_{k+},t)={\sigma }_{xz}^S(x,z_{k-},t)\\ {\sigma }_{}^F(x,z_{k+},t)={\sigma }_{}^F(x,z_{k-},t);{\sigma }_{zz}^I(x,z_{k+},t)={\sigma }_{zz}^I(x,z_{k-},t)\\ {\sigma }_{xz}^I(x,z_{k+},t)={\sigma }_{xz}^I(x,z_{k-},t);u_z^S(x,z_{k+},t)=u_z^S(x,z_{k-},t)\\ u_x^S(x,z_{k+},t)=u_x^S(x,z_{k-},t);u_z^F(x,z_{k+},t)=u_z^F(x,z_{k-},t)\\ u_z^I(x,z_{k+},t)=u_z^I(x,z_{k-},t);u_x^I(x,z_{k+},t)=u_x^I(x,z_{k-},t)\end{array}$$ (5)

In this paper, we consider that when the surface of a layered saturated frozen soil foundation is permeable (z=0), there:

$$\begin{array}{l}{\sigma }_{zz}^S(x\mathrm{,0,}t)+{\sigma }^F(x\mathrm{,0,}t)+{\sigma }_{zz}^I(x\mathrm{,0,}t)=-F(x,t),{\sigma }_{xz}^S(x\mathrm{,0,}t)=\mathrm{0,}\\ {\sigma }^F(x\mathrm{,0,}t)=\mathrm{0,}{\sigma }_{zz}^I(x\mathrm{,0,}t)=\mathrm{0,}{\sigma }_{xz}^I(x\mathrm{,0,}t)=0\end{array}$$ (6)

The boundary conditions are expressed as follows at the bottom impermeable bedrock (z=H):

$$\begin{array}{l}{u}_z^S(x,H,t)=\mathrm{0,}{u}_x^S(x,H,t)=\mathrm{0,}{u}_z^F(x,H,t)=\mathrm{0,}\\ u_x^I(x,H,t)=\mathrm{0,}{u}_z^I(x,H,t)=0\end{array}$$ (7)

2.3 The relationship between temperature and pore ice content

Carcione and Seriani (2001) used the normal distribution pore for the expression of ice content, and they postulated that at a specific temperature when the pore size is below a certain threshold, the water within the pore will not freeze, whereas pores larger than this size will experience complete freezing. Consequently, the correlation between temperature and unfrozen water content is computed as follows:

$${\varphi }_F=\left(1-{\varphi }_S\right)\frac{erf\left(\zeta \right)+erf\left(\eta \right)}{1+erf\left(\eta \right)}$$ (8)

$$\zeta =\frac{r_0/\ln \left(T_0/T_{k0}\right)}{\sqrt{2}\Delta r}-\eta$$ (9)

$$\eta =\frac{r_{av}}{\sqrt{2}\Delta r}$$ (10)

where $r_{av}$ is the average pore radius; $T_{k0}=T+T_0$, $T_{k0}$ and $T$ are temperature, $T_{k0}$ is in Kelvin, $T$ is in Celsius, $T_{k0}$ = 273.15K; The specific values of $r_{av}$, $\Delta r$ and $r_0$ are detailed in Reference (Zhou, 2020).

The volume fraction of each phase in saturated frozen soil as (Zhou and Lai, 2011):

$${\varphi }_S=1-\varphi ,{\varphi }_I=\varphi (1-S_r)$$ (11)

where $S_r$ denotes the saturation level of pore liquid water, $S_r={\varphi }_F/\varphi$; $\varphi$ denotes the porosity of porous media; ${\varphi }_a$ (a=S, F, I) denotes the porosity of a phase.

3 SOLUTION OF THE SATURATED FROZEN SOIL LAYERED SYSTEM

3.1 Solution in a homogeneous layer

The potential functions ${\phi }_a$ and ${\psi }_a\left(a=S,F,I\right)$ are introduced. The displacement can be expressed by Helmholtz vector decomposition theorem as:

$$u^a=\nabla {\phi }_a+\nabla \times {\psi }_a$$ (12)

Substituting Eq. (12) into Eq. (4a), Eq. (4b) and Eq. (4c):

$${\rho }_{11}{\ddot{\phi }}_S+{\rho }_{12}{\ddot{\phi }}_F+{\rho }_{13}{\ddot{\phi }}_I=R_{11}{\nabla }^2{\phi }_S+R_{12}{\nabla }^2{\phi }_F+R_{13}{\nabla }^2{\phi }_I-\left(b_{12}+b_{13}\right){\dot{\phi }}_S+b_{12}{\dot{\phi }}_F+b_{13}{\dot{\phi }}_I$$ (13a)

$${\rho }_{12}{\ddot{\phi }}_S+{\rho }_{22}{\ddot{\phi }}_F+{\rho }_{23}{\ddot{\phi }}_I=R_{12}{\nabla }^2{\phi }_S+R_{22}{\nabla }^2{\phi }_F+R_{23}{\nabla }^2{\phi }_I+b_{12}{\dot{\phi }}_S-\left(b_{12}+b_{23}\right){\dot{\phi }}_F+b_{23}{\dot{\phi }}_I$$ (13b)

$${\rho }_{13}{\ddot{\phi }}_S+{\rho }_{23}{\ddot{\phi }}_F+{\rho }_{33}{\ddot{\phi }}_I=R_{13}{\nabla }^2{\phi }_S+R_{23}{\nabla }^2{\phi }_F+R_{33}{\nabla }^2{\phi }_I+b_{13}{\dot{\phi }}_S+b_{23}{\dot{\phi }}_F-\left(b_{13}+b_{23}\right){\dot{\phi }}_I$$ (13c)

$${\rho }_{11}{\ddot{\psi }}_S+{\rho }_{12}{\ddot{\psi }}_F+{\rho }_{13}{\ddot{\psi }}_I={\mu }_{11}{\nabla }^2{\psi }_S+{\mu }_{13}{\nabla }^2{\psi }_I-\left(b_{12}+b_{13}\right){\dot{\psi }}_S+b_{12}{\dot{\psi }}_F+b_{13}{\dot{\psi }}_I$$ (14a)

$${\rho }_{12}{\ddot{\psi }}_S+{\rho }_{22}{\ddot{\psi }}_F+{\rho }_{23}{\ddot{\psi }}_I=b_{12}{\dot{\psi }}_S-\left(b_{12}+b_{23}\right){\dot{\psi }}_F+b_{23}{\dot{\psi }}_I$$ (14b)

$${\rho }_{13}{\ddot{\psi }}_S+{\rho }_{23}{\ddot{\psi }}_F+{\rho }_{33}{\ddot{\psi }}_I={\mu }_{13}{\nabla }^2{\psi }_S+{\mu }_{33}{\nabla }^2{\psi }_I+b_{13}{\dot{\psi }}_S+b_{23}{\dot{\psi }}_F-\left(b_{13}+b_{23}\right){\dot{\psi }}_I$$ (14c)

Under the action of moving harmonic load, all variables can be expressed as:

$$G(x,z,t)=G^{*}(x,z){\text{e}}^{\text{i}\omega t}$$ (15)

where i is the imaginary unit, for the convenience of representation, the superscript ' * ' in the following text is omitted.

As shown in Figure 1, the coordinate is transformed into $x=x_1-ct$, $z=x_3$ and the potential functions $\phi$ and $\psi$ in the moving coordinate system can be expressed as follows:

$$\chi \left(x_1-ct,x_3,t\right)=\chi \left(x,z\right)e^{i\omega t}$$ (16a)

$$\dot{\chi }\left(x_1-ct,x_3,t\right)=\left(i\omega \chi -c{\chi }_{,x}\right)e^{i\omega t}$$ (16b)

$$\ddot{\chi }\left(x_1-ct,x_3,t\right)=\left(-{\omega }^2\chi -2i\omega c{\chi }_{,x}+c{\chi }_{,xx}\right)e^{i\omega t}$$ (16c)

where $\chi =\phi ,\psi$, ${\chi }_{,x}$ and ${\chi }_{,xx}$ denote the partial differential of $x$.

The Fourier transform of $x$ is:

$$\tilde{f}(\xi ,z)={\displaystyle {\int }_{-\infty }^{\infty }f(x,z)e^{-i\xi x}dx}$$ (17)

where the superscript ' ~ ' indicates that the spatial coordinate x is subjected to a Fourier transform.

Combining Eq. (16a-16c) and Fourier transform of Eq. (13a), Eq. (13b), and Eq. (13c), we can obtain:

$$R_{11}\frac{d^2{\tilde{\phi }}_S}{d{z}^2}+A_{11}{\tilde{\phi }}_S+R_{12}\frac{d^2{\tilde{\phi }}_F}{d{z}^2}+A_{12}{\tilde{\phi }}_F+R_{13}\frac{d^2{\tilde{\phi }}_I}{d{z}^2}+A_{13}{\tilde{\phi }}_I=0$$ (18a)

$$R_{12}\frac{d^2{\tilde{\phi }}_S}{d{z}^2}+A_{12}{\tilde{\phi }}_S+R_{22}\frac{d^2{\tilde{\phi }}_F}{d{z}^2}+A_{22}{\tilde{\phi }}_F+R_{23}\frac{d^2{\tilde{\phi }}_I}{d{z}^2}+A_{23}{\tilde{\phi }}_I=0$$ (18b)

$$R_{13}\frac{d^2{\tilde{\phi }}_S}{d{z}^2}+A_{13}{\tilde{\phi }}_S+R_{23}\frac{d^2{\tilde{\phi }}_F}{d{z}^2}+A_{23}{\tilde{\phi }}_F+R_{33}\frac{d^2{\tilde{\phi }}_I}{d{z}^2}+A_{33}{\tilde{\phi }}_I=0$$ (18c)

where $A_{11}=-R_{11}{\xi }^2-{\rho }_{11}{\alpha }_1-\left(b_{12}+b_{13}\right){\alpha }_2$, $A_{12}=-R_{12}{\xi }^2-{\rho }_{12}{\alpha }_1+b_{12}{\alpha }_2$, $A_{13}=-R_{13}{\xi }^2-{\rho }_{13}{\alpha }_1+b_{13}{\alpha }_2$, $A_{22}=-R_{22}{\xi }^2-{\rho }_{22}{\alpha }_1-\left(b_{12}+b_{23}\right){\alpha }_2$, $A_{23}=-R_{23}{\xi }^2-{\rho }_{23}{\alpha }_1+b_{23}{\alpha }_2$, ${\alpha }_1=-{\omega }^2+2\omega c\xi -c^2{\xi }^2$ , ${\alpha }_2=i\omega -ic\xi$.

Let's assume that the solution of Eqs. (18) can be represented by the following form:

$${\left[{\tilde{\phi }}_S{\tilde{\phi }}_F{\tilde{\phi }}_I\right]}^T={[c^S{c}^F{c}^I]}^T\exp (\lambda z)$$ (19)

By substituting Eq. (19) into Eq. (18a), Eq. (18b) and Eq. (18c), the linear equations are obtained:

$$\left[\begin{array}{ccc}{\lambda }^2{R}_{13}+A_{13}& {\lambda }^2{R}_{12}+A_{12}& {\lambda }^2{R}_{11}+A_{11}\\ {\lambda }^2{R}_{23}+A_{23}& {\lambda }^2{R}_{22}+A_{22}& {\lambda }^2{R}_{12}+A_{12}\\ {\lambda }^2{R}_{33}+A_{33}& {\lambda }^2{R}_{23}+A_{23}& {\lambda }^2{R}_{13}+A_{13}\end{array}\right]\left[\begin{array}{c}{c}^I\\ c^F\\ c^S\end{array}\right]=0$$ (20)

The following equation can be obtained:

$${\beta }_1{\lambda }^6+{\beta }_2{\lambda }^4+{\beta }_3{\lambda }^2+{\beta }_4=0$$ (21)

where ${\beta }_1$, ${\beta }_2$, ${\beta }_3$ and ${\beta }_4$ are shown in ^{Appendix 2}.

Suppose that the root of Eq. (21) is $\pm {\lambda }_n\left(n=\mathrm{1,}\mathrm{2,}3\right)$, then ${\lambda }_n$ are given by:

$${\lambda }_n=\sqrt{d_n}\left(\mathrm{Re}\left[{\lambda }_n\right]\ge \mathrm{0,}n=\mathrm{1,}\mathrm{2,}3\right)$$ (22)

The generalized solution of the potential function can be found as:

$${\tilde{\phi }}_S={\displaystyle \sum _{n=1}^3\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)}$$ (23a)

$${\tilde{\phi }}_F={\displaystyle \sum _{n=1}^3{\delta }_{pn}^F\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)}$$ (23b)

$${\tilde{\phi }}_I={\displaystyle \sum _{n=1}^3{\delta }_{pn}^I\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)}$$ (23c)

where ${\delta }_{pn}^F$, ${\delta }_{pn}^I\left(n=\mathrm{1,}\mathrm{2,}3\right)$ are given in ^{Appendix 2}; $D_1$, $D_2$, $D_3$, $E_1$, $E_2$, and $E_3$ are undetermined coefficients.

The Fourier transform of Eqs. (14a), (14b) and (14c) yields the following:

$${\mu }_{11}\frac{d^2{\tilde{\psi }}_S}{d{z}^2}+B_{11}{\tilde{\psi }}_S+B_{12}{\tilde{\psi }}_F+{\mu }_{13}\frac{d^2{\tilde{\psi }}_I}{d{z}^2}+B_{13}{\tilde{\psi }}_I=0$$ (24a)

$$B_{21}{\tilde{\psi }}_S+B_{22}{\tilde{\psi }}_F+B_{23}{\tilde{\psi }}_I=0$$ (24b)

$${\mu }_{13}\frac{d^2{\tilde{\psi }}_S}{d{z}^2}+B_{31}{\tilde{\psi }}_S+B_{32}{\tilde{\psi }}_F+{\mu }_{33}\frac{d^2{\tilde{\psi }}_I}{d{z}^2}+B_{33}{\tilde{\psi }}_I=0$$ (24c)

where $B_{11}=-{\mu }_{11}{\xi }^2-{\rho }_{11}{\alpha }_1-(b_{12}+b_{13}){\alpha }_2$, $B_{12}=-{\rho }_{12}{\alpha }_1+b_{12}{\alpha }_2$, $B_{13}=-{\mu }_{13}{\xi }^2-{\rho }_{13}{\alpha }_1+b_{13}{\alpha }_2$, $B_{21}={\rho }_{12}{\alpha }_1-b_{12}{\alpha }_2$, $B_{22}={\rho }_{22}{\alpha }_1+\left(b_{12}+b_{23}\right){\alpha }_2$, $B_{23}={\rho }_{23}{\alpha }_1-b_{23}{\alpha }_2$, $B_{31}=-{\mu }_{13}{\xi }^2-{\rho }_{13}{\alpha }_1+b_{13}{\alpha }_2$, $B_{32}=-{\rho }_{23}{\alpha }_1+b_{23}{\alpha }_2$,$B_{33}=-{\mu }_{33}{\xi }^2-{\rho }_{33}{\alpha }_1-\left(b_{13}+b_{23}\right){\alpha }_2$.

Suppose the solution of Eqs. (24) has the following form:

$${\left[{\tilde{\psi }}_S{\tilde{\psi }}_F{\tilde{\psi }}_I\right]}^T={[h^S{h}^F{h}^I]}^T\exp (rz)$$ (25)

By substituting Eq. (25) into Eq. (24a), Eq. (24b) and Eq. (24c), the linear equations are obtained:

$$\left[\begin{array}{ccc}{\mu }_{13}{r}^2+B_{13}& B_{12}& {\mu }_{11}{r}^2+B_{11}\\ B_{23}& B_{22}& B_{21}\\ {\mu }_{33}{r}^2+B_{33}& B_{32}& {\mu }_{13}{r}^2+B_{31}\end{array}\right]\left[\begin{array}{c}{h}^I\\ h^F\\ h^S\end{array}\right]=0$$ (26)

The following equation can be obtained:

$${\beta }_5{r}^4+{\beta }_6{r}^2+{\beta }_7=0$$ (27)

where ${\beta }_5$, ${\beta }_6$ and ${\beta }_7$ are shown in ^{Appendix 2}.

Suppose that the root of Eq. (24) is $\pm r_n\left(n=\mathrm{1,}2\right)$, then $r_n$ are given by:

$$r_n=\sqrt{t_n}\left(\mathrm{Re}\left[r_n\right]\ge \mathrm{0,}n=1,2\right)$$ (28)

The generalized solution of the potential function can be found as:

$${\tilde{\psi }}_S={\displaystyle \sum _{n=1}^2\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (29a)

$${\tilde{\psi }}_F={\displaystyle \sum _{n=1}^2{\delta }_{sn}^F\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (29b)

$${\tilde{\psi }}_I={\displaystyle \sum _{n=1}^2{\delta }_{sn}^I\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (29c)

where ${\delta }_{sn}^F$, ${\delta }_{sn}^I\left(n=\mathrm{1,}\mathrm{2,}3\right)$ are given in ^{Appendix 2}; $M_1$, $M_2$, $N_2$, $N_2$ are undetermined coefficients.

In the right-angle coordinate system, each displacement component can be expressed as the potential function:

$$u_x=\frac{\partial \phi }{\partial x}-\frac{\partial \psi }{\partial z},u_z=\frac{\partial \phi }{\partial z}+\frac{\partial \psi }{\partial x}$$ (30)

Combining Eq. (2) and Eq. (30) and substituting into Eq. (1), the constitutive equation for saturated frozen soil is:

$${\sigma }_{zz}^S=\left(K_1-2{\mu }_{11}/3\right){\nabla }^2{\phi }_S+C_{12}{\nabla }^2{\phi }_F+\left(C_{13}-{\mu }_{13}/3\right){\nabla }^2{\phi }_I+2{\mu }_{11}\left(\frac{{\partial }^2{\phi }_S}{\partial z^2}+\frac{{\partial }^2{\psi }_S}{\partial x\partial z}\right)+{\mu }_{13}\left(\frac{{\partial }^2{\phi }_I}{\partial z^2}+\frac{{\partial }^2{\psi }_I}{\partial x\partial z}\right)$$ (31a)

$${\sigma }_{xz}^S={\mu }_{11}\left(2\frac{{\partial }^2{\phi }_S}{\partial x\partial z}+\frac{{\partial }^2{\psi }_S}{\partial x^2}-\frac{{\partial }^2{\psi }_S}{\partial z^2}\right)+\frac{1}{2}{\mu }_{13}\left(2\frac{{\partial }^2{\phi }_I}{\partial x\partial z}+\frac{{\partial }^2{\psi }_I}{\partial x^2}-\frac{{\partial }^2{\psi }_I}{\partial z^2}\right)$$ (31b)

$${\sigma }^F=C_{12}{\nabla }^2{\phi }_S+K_2{\nabla }^2{\phi }_F+C_{23}{\nabla }^2{\phi }_I$$ (31c)

$${\sigma }_{zz}^I=\left(C_{13}-{\mu }_{13}/3\right){\nabla }^2{\phi }_S+C_{23}{\nabla }^2{\phi }_F+\left(K_3-2{\mu }_{33}/3\right){\nabla }^2{\phi }_I+{\mu }_{13}\left(\frac{{\partial }^2{\phi }_S}{\partial z^2}+\frac{{\partial }^2{\psi }_S}{\partial x\partial z}\right)+2{\mu }_{33}\left(\frac{{\partial }^2{\phi }_I}{\partial z^2}+\frac{{\partial }^2{\psi }_I}{\partial x\partial z}\right)$$ (31d)

$${\sigma }_{xz}^I={\mu }_{33}\left(2\frac{{\partial }^2{\phi }_I}{\partial x\partial z}+\frac{{\partial }^2{\psi }_I}{\partial x^2}-\frac{{\partial }^2{\psi }_I}{\partial z^2}\right)+\frac{1}{2}{\mu }_{13}\left(2\frac{{\partial }^2{\phi }_S}{\partial x\partial z}+\frac{{\partial }^2{\psi }_S}{\partial x^2}-\frac{{\partial }^2{\psi }_S}{\partial z^2}\right)$$ (31d)

Fourier transform is applied to Eq. (30) and Eq. (31):

$${\tilde{u}}_x=\text{i}\xi \tilde{\phi }-\frac{\partial \tilde{\psi }}{\partial z},{\tilde{u}}_z=\frac{\partial \phi }{\partial z}+\text{i}\xi \tilde{\psi }$$ (32)

$$\begin{array}{l}{\tilde{\sigma }}_{zz}^S=\left(K_1-2{\mu }_{11}/3\right)\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_S+C_{12}\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_F+\left(C_{13}-{\mu }_{13}/3\right)\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_I\\ +2{\mu }_{11}\left(\frac{{\partial }^2{\tilde{\phi }}_S}{\partial z^2}+\text{i}\xi \frac{\partial {\tilde{\psi }}_S}{\partial z}\right)+{\mu }_{13}\left(\frac{{\partial }^2{\tilde{\phi }}_I}{\partial z^2}+\text{i}\xi \frac{\partial {\tilde{\psi }}_I}{\partial z}\right)\end{array}$$ (33a)

$${\tilde{\sigma }}_{xz}^S={\mu }_{11}\left(2\text{i}\xi \frac{\partial {\tilde{\phi }}_S}{\partial z}-{\xi }^2{\tilde{\psi }}_S-\frac{{\partial }^2{\tilde{\psi }}_S}{\partial z^2}\right)+\frac{1}{2}{\mu }_{13}\left(2\text{i}\xi \frac{\partial {\phi }_I}{\partial z}-{\xi }^2{\tilde{\psi }}_I-\frac{{\partial }^2{\tilde{\psi }}_I}{\partial z^2}\right)$$ (33b)

$${\tilde{\sigma }}^F=C_{12}\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_S+K_2\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_F+C_{23}\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_I$$ (33c)

$$\begin{array}{l}{\tilde{\sigma }}_{zz}^I=\left(C_{13}-{\mu }_{13}/3\right)\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_S+C_{23}\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_F+\left(K_3-2{\mu }_{33}/3\right)\left(-{\xi }^2+\frac{{\partial }^2}{\partial z^2}\right){\tilde{\phi }}_I+\\ {\mu }_{13}\left(\frac{{\partial }^2{\tilde{\phi }}_S}{\partial z^2}+\text{i}\xi \frac{\partial {\tilde{\psi }}_S}{\partial z}\right)+2{\mu }_{33}\left(\frac{{\partial }^2{\tilde{\phi }}_I}{\partial z^2}+\text{i}\xi \frac{\partial {\tilde{\psi }}_I}{\partial z}\right)\end{array}$$ (33d)

$${\tilde{\sigma }}_{xz}^I={\mu }_{33}\left(2\text{i}\xi \frac{\partial {\tilde{\phi }}_I}{\partial z}-{\xi }^2{\tilde{\psi }}_I-\frac{{\partial }^2{\tilde{\psi }}_I}{\partial z^2}\right)+\frac{1}{2}{\mu }_{13}\left(2\text{i}\xi \frac{\partial {\tilde{\phi }}_S}{\partial z}-{\xi }^2{\tilde{\psi }}_S-\frac{{\partial }^2{\tilde{\psi }}_S}{\partial z^2}\right)$$ (33e)

By substituting Eqs. (23) and (29) into Eqs. (32)-(33), expressions for the displacements and stresses of the phases in the saturated frozen soil medium in the Fourier transform domain are obtained:

$${\tilde{u}}_x^S=i\xi {\displaystyle \sum _{n=1}^3\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)}+{\displaystyle \sum _{n=1}^2{r}_n\left(M_n{e}^{-r_{n}z}-N_n{e}^{r_{n}z}\right)}$$ (34a)

$${\tilde{u}}_z^S={\displaystyle \sum _{n=1}^3-{\lambda }_n\left(D_n{e}^{-{\lambda }_{n}z}-E_n{e}^{{\lambda }_{n}z}\right)}+i\xi {\displaystyle \sum _{n=1}^2\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (34b)

$${\tilde{u}}_z^F={\displaystyle \sum _{n=1}^3-{\lambda }_n{\delta }_{pn}^F\left(D_n{e}^{-{\lambda }_{n}z}-E_n{e}^{{\lambda }_{n}z}\right)}+i\xi {\displaystyle \sum _{n=1}^2{\delta }_{sn}^F\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (34c)

$${\tilde{u}}_x^I=i\xi {\displaystyle \sum _{n=1}^3{\delta }_{pn}^I\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)}+{\displaystyle \sum _{n=1}^2{r}_n{\delta }_{sn}^I\left(M_n{e}^{-r_{n}z}-N_n{e}^{r_{n}z}\right)}$$ (34d)

$${\tilde{u}}_z^I={\displaystyle \sum _{n=1}^3-{\lambda }_n{\delta }_{pn}^I\left(D_n{e}^{-{\lambda }_{n}z}-E_n{e}^{{\lambda }_{n}z}\right)}+i\xi {\displaystyle \sum _{n=1}^2{\delta }_{sn}^I\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (34e)

$${\tilde{\sigma }}_{xz}^S=i\xi {\displaystyle \sum _{n=1}^3-{\lambda }_n\left(2{\mu }_{11}+{\mu }_{13}{\delta }_{pn}^I\right)\left(D_n{e}^{-{\lambda }_{n}z}-E_n{e}^{{\lambda }_{n}z}\right)}-{\displaystyle \sum _{n=1}^2\left(r_n{}^2+{\xi }^2\right)\left({\mu }_{11}+\frac{1}{2}{\mu }_{13}{\delta }_{sn}^I\right)\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (34f)

$${\tilde{\sigma }}_{zz}^S={\displaystyle \sum _{n=1}^3{f}_{Sn}\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)}-i\xi {\displaystyle \sum _{n=1}^2{r}_n\left(2{\mu }_{11}+{\mu }_{13}{\delta }_{sn}^I\right)\left(M_n{e}^{-r_{n}z}-N_n{e}^{r_{n}z}\right)}$$ (34g)

$${\tilde{\sigma }}_{}^F={\displaystyle \sum _{n=1}^3\left(R_{12}+R_{22}{\delta }_{pn}^F+R_{23}{\delta }_{pn}^I\right)}\left({\lambda }_n^2-{\xi }^2\right)\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)$$ (34h)

$${\tilde{\sigma }}_{xz}^I=i\xi {\displaystyle \sum _{n=1}^3-{\lambda }_n\left(2{\mu }_{33}{\delta }_{pn}^I+{\mu }_{13}\right)\left(D_n{e}^{-{\lambda }_{n}z}-E_n{e}^{{\lambda }_{n}z}\right)}-{\displaystyle \sum _{n=1}^2\left(r_n{}^2+{\xi }^2\right)\left({\mu }_{33}{\delta }_{sn}^I+\frac{1}{2}{\mu }_{13}\right)\left(M_n{e}^{-r_{n}z}+N_n{e}^{r_{n}z}\right)}$$ (34i)

$${\tilde{\sigma }}_{zz}^I={\displaystyle \sum _{n=1}^3{f}_{In}\left(D_n{e}^{-{\lambda }_{n}z}+E_n{e}^{{\lambda }_{n}z}\right)}-i\xi {\displaystyle \sum _{n=1}^2{r}_n\left({\mu }_{13}+2{\mu }_{33}{\delta }_{sn}^I\right)\left(M_n{e}^{-r_{n}z}-N_n{e}^{r_{n}z}\right)}$$ (34j)

where $f_{Sn}=\left[\left(R_{11}-2{\mu }_{11}\right)+R_{12}{\delta }_{pn}^F+\left(R_{13}-{\mu }_{13}\right){\delta }_{pn}^I\right]\left({\lambda }_n^2-{\xi }^2\right)+\left(2{\mu }_{11}+{\mu }_{13}{\delta }_{pn}^I\right){\lambda }_n^2\left(n=\mathrm{1,}\mathrm{2,}3\right)$, $f_{In}=\left[\left(R_{13}-{\mu }_{13}\right)+R_{23}{\delta }_{pn}^F+\left(R_{33}-2{\mu }_{33}\right){\delta }_{pn}^I\right]\left({\lambda }_n^2-{\xi }^2\right)+\left({\mu }_{13}+2{\mu }_{33}{\delta }_{pn}^I\right){\lambda }_n^2\left(n=\mathrm{1,}\mathrm{2,}3\right)$.

3.2 Solution of the saturated frozen soil layered system

Based on the analysis conducted in the previous section, the stress and displacement in the $k$ layer of saturated frozen soil can be expressed in the following matrix form:

$$\left[\begin{array}{c}{L}^k\left(\xi ,z\right)\\ N^k\left(\xi ,z\right)\end{array}\right]=\left[S^k\right]E^k\left(z\right)\left[R^k\right]$$ (35)

where the coefficient superscript '$k$' represents the $k$ layer of the soil, $L\left(z\right)={\left[\begin{array}{cccc}{\tilde{\sigma }}_{zz}^S+{\tilde{\sigma }}_{}^F+{\tilde{\sigma }}_{xz}^I& {\tilde{\sigma }}_{xz}^S& {\tilde{\sigma }}_{}^F& {\tilde{\sigma }}_{zz}^I\end{array}{\tilde{\sigma }}_{xz}^I\right]}^{\text{T}}$, $N\left(z\right)={\left[\begin{array}{cccc}{\tilde{u}}_x^S& {\tilde{u}}_z^S& {\tilde{u}}_z^F& {\tilde{u}}_z^I\end{array}{\tilde{u}}_z^I\right]}^{\text{T}}$;${\left[S\right]}_{10\times 10}$ is shown in ^{Appendix 3}, $\left[R\right]=[\begin{array}{cccc}{D}_1& D_2& D_3& M_1\end{array}\begin{array}{cccccc}{M}_2& E_1& E_2& E_3& N_1& N_2{]}^T\end{array}$, $E\left(z\right)=\text{diag[}\begin{array}{cccccc}{e}^{-{\lambda }_{1}z}& e^{-{\lambda }_{2}z}& e^{-{\lambda }_{3}z}& e^{-r_{1}z}& e^{-r_{2}z}& e^{{\lambda }_{1}z}\end{array}\begin{array}{cccc}{e}^{{\lambda }_{2}z}& e^{{\lambda }_{3}z}& e^{r_{1}z}& e^{r_{2}z}\end{array}]$, $D_1$, $D_2$, $D_3$, $M_1$, $M_2$, $E_1$, $E_2$, $E_3$, $N_1$, $N_2$ are undetermined coefficients.

Then the stress and displacement of the adjacent $k$+1 layers can be written as similar to Eq.(35):

$$\left[\begin{array}{c}{L}^{k+1}\left(\xi ,z\right)\\ N^{k+1}\left(\xi ,z\right)\end{array}\right]=\left[S^{k+1}\right]E^{k+1}\left(z\right)\left[R^{k+1}\right]$$ (36)

To establish a local coordinate system, the following continuity conditions must be satisfied between adjacent soil layers in saturated frozen soils:

$${\left.\left[\begin{array}{c}{L}^{k+1}\left(\xi ,z\right)\\ N^{k+1}\left(\xi ,z\right)\end{array}\right]\right|}_{z=h^k}={\left.\left[\begin{array}{c}{L}^k\left(\xi ,z\right)\\ N^k\left(\xi ,z\right)\end{array}\right]\right|}_{z=0}$$ (37)

where $h^k=z_{k+1}-z_k$.

Substituting Eq. (36) into Eq. (37) gives:

$$\left[R^{k+1}\right]={\left[S^{k+1}\right]}^{-1}\left[S^k\right]E^k\left(h^k\right)\left[R^k\right]=\left[T^k\right]\left[R^k\right]$$ (38)

where ${\left[T^k\right]}_{10\times 10}$ is the transfer matrix.

According to the recurrence Eq. (38), the amplitude vector of any layer in saturated frozen soil can be related to the amplitude vector of the first layer:

$$\left[R^{k+1}\right]=\left[T^k\right]\left[T^{k-1}\right]\cdots \left[T^1\right]\left[R^1\right]$$ (39)

The time domain expression for the moving harmonic load is given below:

$$F\left(x,t\right)=\left\{\begin{array}{l}q{\text{e}}^{\text{i}\omega t}=\frac{q_0}{2l}{\text{e}}^{\text{i}\omega t}-l<x<l\\ 0\text{else}\end{array}\right.$$ (40)

Fourier transform of the load:

$$F\left(\xi ,t\right)={\displaystyle {\int }_{-l}^l\frac{q_0}{2l}{\text{e}}^{\text{i}\omega t}{\text{e}}^{-j\xi x}}dx=\frac{q_0{\text{e}}^{\text{i}\omega t}}{2l}{\displaystyle {\int }_{-l}^l{e}^{-j\xi x}dx=}\frac{q_0{\text{e}}^{\text{i}\omega t}}{2l}{\displaystyle {\int }_{-l}^l\left(\cos \xi x-j\sin \xi x\right)dx=}\frac{q_0{\text{e}}^{\text{i}\omega t}\sin \left(\xi l\right)}{\xi l}$$ (41)

Combined Eq. (15), for moving simple harmonic loads subjected to a saturated frozen soil foundation surface with a rigid underlayment, the boundary conditions in the frequency domain can be obtained as:

At z=0:

$$\begin{array}{l}{\tilde{\sigma }}_{zz}^S(\xi \mathrm{,0})+{\tilde{\sigma }}^F(\xi \mathrm{,0})+{\tilde{\sigma }}_{zz}^I(\xi \mathrm{,0})=-q_0\frac{\sin \left(\xi l\right)}{\xi l},{\tilde{\sigma }}_{xz}^S(\xi \mathrm{,0})=\mathrm{0,}\\ {\tilde{\sigma }}^F(\xi \mathrm{,0})=\mathrm{0,}{\tilde{\sigma }}_{zz}^I(\xi \mathrm{,0})=\mathrm{0,}{\tilde{\sigma }}_{zz}^I(\xi \mathrm{,0})=0\end{array}$$ (42)

At z=H:

$${\tilde{u}}_z^S(\xi ,H)=\mathrm{0,}{\tilde{u}}_z^F(\xi ,H)=\mathrm{0,}{\tilde{u}}_z^I(\xi ,H)=\mathrm{0,}{\tilde{u}}_x^I(\xi ,H)=0�{\tilde{u}}_x^S(\xi ,H)=0$$ (43)

Substitute the formula (34) into the boundary condition formula (42) and the formula (43), and the arrangement can be obtained:

$$\left[M^1\right]\left[R^1\right]={\left[\begin{array}{cccc}-q_0\frac{\sin \left(\xi l\right)}{\xi l}& 0& 0& 0\end{array}0\right]}^{\text{T}}$$ (44a)

$$\left[M^N\right]E^N\left(h_N\right)\left[R^N\right]={\left[\begin{array}{cccc}0& 0& 0& 0\end{array}0\right]}^{\text{T}}$$ (44b)

where ${\left[M^1\right]}_{5\times 10}$ and ${\left[M^N\right]}_{5\times 10}$ are detailed in ^{Appendix 4}.

The simultaneous Eq. (39) and Eq. (44a) can be obtained as follows:

$$\left[Q\right]\left[R^j\right]={\left[\begin{array}{cccc}-q_0\frac{\sin \left(\xi l\right)}{\xi l}& 0& 0& 0\end{array}0\right]}^{\text{T}}$$ (45)

where $\left[Q\right]=\left[M^1\right]{\left[T^1\right]}^{-1}\cdots {\left[T^{j-1}\right]}^{-1}$.

The simultaneous Eq. (44b) and Eq. (45) can be obtained as follows:

$${\left[Q{M}^N{E}^N\left(h_N\right)\right]}^T\left[R^N\right]={\left[\begin{array}{cccc}-q_0\frac{\sin \left(\xi l\right)}{\xi l}& 0& 0& 0\end{array}000000\right]}^{\text{T}}$$ (46)

By solving the system of Eq. (46), ten unknown coefficients can be obtained. Combined with Eq. (34), the stress and displacement of any point in the frequency domain can be obtained. Subsequently, the inverse Fourier transform is employed to obtain the corresponding values of stress and displacement in the time domain.

4 NUMERICAL RESULTS AND DISCUSSION

4.1 Verification

Tabatabaie et al. (1994) obtained the analytical solution of the vertical displacement of the soil when the harmonic load is applied to the surface of the single-layer saturated soil foundation by Fourier transform and compared it with the results of the finite element solution, which is in good agreement. In order to verify the accuracy of the numerical calculation in this paper, the dynamic response of layered saturated frozen soil under moving load can be reduced to the dynamic response of single-layer homogeneous saturated soil foundation under harmonic load. The speed of moving harmonic load is taken as 0m/s, and the same soil parameters as saturated soil are set in the layered soil model of saturated frozen soil. The volume fraction of pore ice in saturated frozen soil is taken as 0.001, and the contact parameter is taken as 1, so as to minimize the influence of pore ice content in saturated frozen soil on the dynamic response, so that the layered saturated frozen soil model is degraded into a single-layer homogeneous saturated soil model. It can be compared with the dynamic response solution of two-dimensional saturated soil under harmonic load obtained by Tabatabaie et al. (1994). In order to facilitate comparison with the literature, take H=15m, H/l=15, ${\rho }_S$=2×10³kg/m³, ${\rho }_F$=1×10³kg/m³, ${\mu }_S$=3.85×10⁷N/m², $\varphi$=0.4, ${\varphi }_I$=0.001, $S_r$=0.99, $v$=0.3, $q_0$=1kpa, $\omega$=1rad/s, $c$=0m/s, $\epsilon$=1. Figure 2 shows the variation of vertical displacement with the horizontal direction and the comparison results with the literature (Tabatabaie et al., 1994). It can be seen that the solution in this paper has a high degree of conformity with the literature solution, which verifies the correctness of the theoretical analysis and calculation method in this paper.

Figure 2
Comparison of changes in vertical displacement along the horizontal with Tabatabaie’s solution.

Taking the two-layer foundation as an example, two typical layered saturated frozen soil foundation models of top-soft and bottom-hard and top-hard and bottom-soft are selected to study the influence of surface soil shear modulus, temperature, load moving speed, contact parameters, and load frequency on the dynamic response of the foundation. In the numerical example, the uniform load amplitude $q_0$=1kpa, the load distribution length l=1m, and the thickness of the first layer of soil and the second layer of soil are h₁= 4m and h₂=8m, respectively. Taking the research method of Liu et al. (2022) on the dynamic response of top-soft and bottom-hard and top-hard and bottom-soft foundation, the shear modulus of the surface soil is changed, and the shear modulus of the second layer of soil remains unchanged. According to the different shear modulus of the surface layer of the soil layer, it is divided into five working conditions, that is, working conditions 1, 2, 3, 4, 5 are ${\mu }_{S1}$=6.85GPa, ${\mu }_{S1}$=13.7GPa, ${\mu }_{S1}$=20.55GPa, ${\mu }_{S1}$=27.4GPa, ${\mu }_{S1}$=34.25GPa, respectively. The shear modulus ${\mu }_{S2}$=20.55GPa of the second layer of soil remains unchanged under the five working conditions. Among them, Cases 1, and 2 belong to the top-soft and bottom-hard foundation, Case 3 belongs to the homogeneous foundation, and Cases 4,5 belong to the top-hard and bottom-soft foundation. The double-layer foundation model in the example is shown in Figure 3. Other physical and mechanical parameters of the upper and lower layers of saturated frozen soil are the same as those in Table 1 (Qiu et al., 2018; Li et al., 2019).

Figure 3
Schematic diagram of double-layer foundation model

Table 1
Physical and mechanical parameters of saturated frozen soil foundation

Material parameters	Hydrodynamic Viscosity coefficient	Density of soil particle	Density of liquid	Density of ice particles	Bulk modulus of soil particles	Bulk modulus of liquid	Bulk modulus of ice particle
Symbol (unit)	${\eta }_F$ Kg/(m·s))	${\rho }_S$ (Kg/m3)	${\rho }_F$ (Kg/m3)	${\rho }_I$ (Kg/m3)	$K_S$ (GPa)	$K_F$ (GPa)	$K_I$ (GPa)
Magnitude	1.8×10-3	2580	1000	920	20.9	2.25	8.58
Material parameters	Shear modulus of soil particle	Shear modulus of ice particle	Reference value of soil skeleton permeability coefficient	Reference value of ice permeability coefficient
Symbol (unit)	${\mu }_S$ (GPa)	${\mu }_I$ (GPa)	${\kappa }_{S0}$ (m2)	${\kappa }_{I0}$ (m2)
Magnitude	6.85-27.4	3.32	1.0×10-11	5.0×10-5

4.2 Influence of soil shear modulus

In order to study the influence of soil shear modulus on the dynamic response of saturated frozen soil foundation subjected to moving load, other parameters are kept unchanged (c=60m/s, T₁=T₂= −0.5°C, ${\epsilon }_1$=${\epsilon }_2$= 0.5, f=10Hz). Figure 4 is the vertical displacement and pore water pressure of saturated frozen soil at z=4m with horizontal direction under different working conditions. It can be seen from Figure 4 that the vertical displacement curve is symmetrical about x=0, and the curve amplitude is directly below the action point. There are two critical points x=-3 and x=3 on both sides of the vertical displacement curve. Between the two critical points (-3< x <3), the vertical displacement decreases with the increase of soil shear modulus, while at both ends of the critical point (x >3 and x<-3), the vertical displacement increases with the increase of soil shear modulus. As the soil shear modulus increases, the amplitude of the pore water stress curve decreases gradually, and the amplitude point of the pore water stress curve gradually moves to the negative half-axis of the x-axis. The amplitude of the vertical displacement and pore water stress curve decreases with the increase of the shear modulus of the soil, which is due to the increase of the shear modulus of the surface soil, that is, the rigidity of the soil skeleton of the saturated permafrost soil is enhanced. On the one hand, the saturated frozen soil becomes more rigid, and the deformation of the soil skeleton decreases, resulting in the reduction of vertical displacement. On the other hand, the rigidity of the soil skeleton increases, the load borne by the soil particle phase in saturated frozen soil increases accordingly, and the load borne by the pore water and pore ice decreases. Therefore, the reinforcement of the surface soil in engineering practice can effectively reduce the vertical displacement of the foundation. In addition, when the top-soft and bottom-hard foundation tends to be uniform, that is, Case 1-Case 3, the dynamic response of saturated permafrost foundation decreases gradually, but when the top-hard and bottom-soft foundation tends to be uniform, that is, Case 3-Case 5, the dynamic response of saturated frozen soil foundation increases gradually, which shows that the arrangement order of soft and hard soil layers exerts a significant influence on the horizontal dynamic response of the foundation.

Figure 4
The influence of shear modulus of topsoil on (a) vertical displacement (b) pore water stress at z=4m

Figure 5 is the curve of vertical displacement and pore water stress at x=0 with depth when changing the shear modulus of surface soil. It can be observed that the maximum value of the vertical displacement curve is at the foundation surface, and the vertical displacement continues to decrease with increasing depth, while the displacement curve then converges at the interface between the soil layers and then decreases synchronously. When the softer soil layer is located at the surface, the vertical displacement decays faster with depth than when the harder soil layer is at the surface, indicating that the stress concentration occurs in the top-soft and bottom-hard foundation, while the stress dispersion occurs in the top-hard and bottom-soft foundation. The pore water stress has a maximum value at z=0.8m, and its trend is to increase sharply from 0 near the surface and then decrease rapidly and finally increase slowly, showing an overall increasing trend, which is the same as the change of pore water pressure studied by Liu et al.(2022). The displacement and pore water stress curves of homogeneous foundation (Case 3) are very different from those of inhomogeneous foundation (Case 1, 2, 4, 5) at the interface of the soil layer (z=4m). The variation in shear modulus of the surface soil significantly affects the vertical dynamic response of the foundation.

Figure 5
The influence of shear modulus of surface soil on (a) vertical displacement (b) pore water stress along the depth direction

4.3 Influence of load moving speed

To facilitate the study of the effect of load movement velocity on the dynamic response of layered saturated frozen soil foundation, the dimensionless velocity N₀=c/c_s is defined by referring to the treatment of velocity by Ma et al ( 2023a ), where c_s is the shear wave velocity of surface soil and c is the load moving speed. When the shear modulus of surface soil is 6.85GPa, c_s=1280m/s can be obtained according to Reference (Ma et al., 2023b). Figure 6 and Figure 7 illustrate the curves of vertical displacement and pore water stress at z=4m with dimensionless velocity N₀=0.5, 0.8, 1.0, and 1.2 under Cases 1 and 5, and other parameters are kept unchanged (c=60m/s, T₁=T₂= −0.5°C, ${\epsilon }_1$=${\epsilon }_2$= 0.5, f=10Hz). As depicted in Figure 6, the displacement and pore water stress exhibit an initial increase followed by a subsequent decrease in the case of top-soft and bottom-hard foundation. At low speeds (i.e., N₀=0.5, 0.8) the increase in displacement and pore water stress magnitude with increasing moving speed is small. When the moving speed closely approaches the shear wave velocity of the surface soil (i.e., N₀=1), there is a sharp increase in both the amplitude of displacement and pore water stress, and the displacement and pore water stress of the soil appears the maximum value, indicating that the soil resonance is easy to occur when the moving speed is close to the shear wave velocity of the soil layer. When the moving speed exceeds the shear wave velocity (i.e., N₀=1.2), the vertical displacement amplitude and pore water stress amplitude decrease with continued increase in velocity, but the vertical displacement amplitude at this point is still larger than the vertical displacement amplitude at low speed, and we observed a sharp decrease in the amplitude of pore water stress, even smaller than that at low speed (N₀=0.8). In contrast, in the top-hard and bottom-soft foundation, the amplitude of vertical displacement continues to increase with the increase of the load movement rate, and the amplitude of its growth also increases. This phenomenon can be attributed to the fact that the surface soil modulus of the top-hard and bottom-soft foundation is larger compared with the top-soft and bottom-hard foundation, and the shear wave velocity of the topsoil increases with the increase of the shear modulus of the topsoil, so the vertical displacement does not change similar to that of the top-soft and bottom-hard foundation. The pore water stress varies with the rate of load moving similarly as in the case of top-soft and bottom-hard foundations. The amplitude of pore water stress increases first and then decreases with the increase in speed. Different from the top-soft and bottom-hard foundation, the maximum value of pore water stress in the top-hard and bottom-soft foundation appears at N₀=0.8. In addition, it can be found that the increase of the moving speed changes the pore water stress in the foundation and the vibration phase of the displacement curve in the top-soft and bottom-hard foundation, making the curve more and more complex.

Figure 6
The influence of load moving speed on (a) vertical displacement (b) pore water stress of foundation at z=4m in Case 1.

Figure 7
The influence of load moving speed on (a) vertical displacement (b) pore water stress of foundation at z=4m in Case 5.

4.4 Influence of temperature

From the Reference (Jamin et al., 2021), it is known that soil temperature in frozen areas is seasonally variable above a certain depth, with the upper soil temperature being higher than the lower layer in summer and the upper soil temperature being lower than the lower layer in winter, and that the temperature change also affects the changes in the components in the saturated frozen soils, so the influence of temperature on the dynamic response can not be ignored. Liu and Zhang (2012) pointed out that the unfrozen water content and pore ice content of frozen soil change very sharply at T=0°C ~ -1°C. Even a small change in temperature will cause a large change in the unfrozen water content, and refer to the value of Li et al. (2019) on the temperature of saturated frozen soil, the temperature range studied in this paper is finally determined to be T=-0.3°C ~ -0.9°C. To investigate the impact of temperature on the dynamic response of both top-soft and bottom-hard foundations and top-hard and bottom-soft foundations, and other parameters are kept unchanged (c=60m/s, ${\epsilon }_1$=${\epsilon }_2$= 0.5, f=10Hz). Figure 8 shows the vertical displacement and pore water stress along the horizontal direction at z=4m under the condition of changing the surface soil temperature in top-soft and bottom-hard foundation (Case 1) and top-hard and bottom-soft foundation (Case 5). It can be seen from the figures that the displacement amplitude and pore water stress amplitude increase with increasing temperature in both working conditions. This is because as the temperature increases, the content of ice particles in the pores gradually decreases, and the interaction with the soil particle skeleton gradually weakens, the load it bears is correspondingly reduced, resulting in an increase in the load borne by the pore water and the soil skeleton, which in turn leads to an increase in the vertical displacement of the soil skeleton. By observing the magnitude of changes in the magnitude of vertical displacement and pore water pressure, it can be found that the increase in the magnitude of vertical displacement under Case 1 is smaller than that of vertical displacement under Case 5, while the increase in the magnitude of pore water pressure in Case 1 is larger than that in Case 5, which indicates that the vertical displacement in the top-hard and bottom-soft foundation and pore water pressure in the top-soft and bottom-hard foundation are more sensitive to changes in temperature. In addition, in the top-hard and bottom-soft foundation, it can be observed that when the temperature of the surface soil approaches the lower soil (Case 5, T₁=-0.3°C), there is only a slight increase in the amplitude of the pore water stress curve, and the point of maximum amplitude in the curve shifts towards the negative half of the x-axis.

Figure 8
The influence of temperature on (a) vertical displacement (b) pore water stress at z=4m under different working conditions.

4.5 Influence of contact parameter

The contact parameter $\epsilon$, ranging from 0 to 1, represents the extent to which the soil particle skeleton in saturated frozen soil is supported by the pore ice. When $\epsilon$=1, the ice is suspended in the pores, and the ice has the smallest support to the soil particle skeleton; when $\epsilon$=0 the ice has the largest support to the soil particles. The contact parameter affects the loads carried by each phase of saturated frozen soil by influencing the modulus of the skeleton. Keep other parameters unchanged (c=60m/s, T₁=T₂= −0.5°C, f=10Hz), Figure 9 analyzes the effect of different contact parameters in the surface layer on the dynamic response at z=4m for different working conditions. The results show that vertical displacement and pore water stress amplitude in the three conditions are increased with the increase of contact parameters, the reason for this phenomenon is because with the increase of contact parameters, pore ice particles and soil skeleton of the contact is gradually reduced, the friction between the pore ice particles and soil particles is reduced, the relative movement amplitude in the gradual increase, so the pore water and soil particles skeleton bear the load increases, which in turn cause the foundation vertical displacement and pore water stress increase. The increase of the vertical displacement amplitude of the top-soft and bottom-hard hard foundation is less than that of the homogeneous foundation and the top-hard and bottom-soft foundation, while the increase of the pore water stress of the upper soft and lower hard foundation is just the opposite. In addition, it can be observed that the change of contact parameters has little effect on the vertical displacement, the possible reason is that the ice content in the saturated frozen soil is small at this temperature (T₁=T₂=-0.5°C), the solid ice does not fill the pores in the saturated frozen soil, and soil particles and pore water carry most of the loads, the pore ice particles bear very small stresses, so changing the contact parameters has little effect on the vertical displacement.

Figure 9
The influence of contact parameters on (a) vertical displacement (b) pore water stress at z=4m under different working conditions.

4.6 Influence of load frequency

The specification (Technical Standard of Highway Engineering, 2015) points out that the frequency range of traffic moving load in road engineering is 0Hz-50Hz, so f=0Hz-50Hz is taken as the research range of load frequency in this paper. Keep other parameters unchanged (c=60m/s, T₁=T₂= −0.5°C, ${\epsilon }_1$=${\epsilon }_2$= 0.5), Figure 10 analyzes the variation of dynamic response with load frequency at moving coordinate system (0,4) for different working conditions. The results indicate that the vertical displacement exhibits noticeable fluctuations with changes in load frequency, and the homogeneous foundation (i.e., Case 3) demonstrates maximum amplitude at approximately 18Hz. This is because the wave generated under vibration propagates and reflects between the foundation and the bedrock so that the vibration amplitude in some areas is strengthened, while the vibration amplitude in some areas is weakened. The vertical displacement of the foundation exhibits an increasing trend with the increase in load frequency. The pore water stress curve changes significantly in 0Hz-20Hz, and the amplitude of the curve in Case 4 reaches the maximum at about 8Hz. As the frequency gradually increases (20Hz-50Hz), the influence of the change of load frequency on the pore water stress in the inhomogeneous foundation gradually weakens, while the pore water stress in the homogeneous foundation is less affected by the change of load frequency. As a whole, the vertical displacement is greater for top-soft and bottom-hard foundations compared to top-hard and bottom-soft foundations, while the pore water stress is smaller than that of the top-hard and bottom-soft foundation, which aligns with the conclusion depicted in Figure 4.

Figure 10
The influence of load frequency on (a) vertical displacement (b) pore water stress in different layered foundation at (0,4).

5 CONCLUSIONS

In this paper, the dynamic response of saturated frozen soil foundations under the influence of moving loads is investigated using the transfer matrix method. The analysis of shear modulus, temperature, load moving speed, contact parameters, and load frequency on the displacement and stress of a layered saturated frozen soil foundation yields the following main conclusions:

1. The arrangement order of the soft and hard soil layers within the foundation soil layer has a notable influence on both vertical displacement and pore water stress. In the horizontal direction, the amplitude of the vertical displacement and pore water stress curves decreases as the shear modulus of the surface soil increases. In the vertical direction, the vertical displacement decreases with the increase of depth, and the pore water stress increases sharply near the surface and then decreases rapidly and then increases slowly. The changes of vertical displacement and pore water stress at the soil interface of the top-soft and bottom-hard foundation and the top-hard and bottom-soft foundation are very different from those of homogeneous foundation.

2. With the increase of the load moving speed, the vertical displacement and pore water stress of the top-soft and bottom-hard foundation increase first and then decrease. When the speed is close to the shear wave velocity of the surface soil (N₀=1), the vertical displacement and pore water stress reach the maximum value. In the top-hard and bottom-soft foundation, the vertical displacement increases with the increase of velocity, and the pore water stress increases first and then decreases with the increase of velocity. When N₀=0.8, the pore water stress has the maximum value.

3. In both the top-soft and bottom-hard foundation, as well as the top-hard and bottom-soft foundation, the amplitude of vertical displacement and pore water stress increases with an increase in the temperature of the surface soil. The vertical displacement in the top-hard and bottom-soft foundation and the pore water stress in the top-soft and bottom-hard foundation are more sensitive to the change of temperature.

4. The amplitude of both vertical displacement and pore water stress increases with the contact parameter, and the increase in vertical displacement amplitude for the top-soft and bottom-hard foundation is smaller than the increase in vertical displacement for the homogeneous foundation and the top-hard and bottom-soft foundation, whereas the increase in pore water stress amplitude for the top-soft and bottom-hard foundation is the opposite.

5. The vertical displacement of the foundation increases as the load frequency increases. The pore water stress of the foundation is significantly influenced by the load frequency in the range of 0Hz-20Hz, but with increasing loading frequency (20Hz-50Hz), the influence gradually weakens. On the whole, the vertical displacement of the top-soft and bottom-hard foundation is larger than that of the top-hard and bottom-soft foundation, while the pore water stress is smaller than that of the top-hard and bottom-soft foundation.

Appendix 1

Coefficients $${\mu }_{\text{11}}$$, $${\mu }_{\text{13}}$$, $${\mu }_{33}$$, $$C_{12}$$, $$C_{13}$$, $$C_{23}$$,$$K_1$$, $$K_2$$, $$K_3$$, $$R_{ij}$$ in Eq. (1), Eq. (2) and Eq.(4)

$$\left.\begin{array}{l}{\mu }_{11}={\left[\left(1-g_1\right){\varphi }_S\right]}^2{\mu }_{av}+{\mu }_{sm}\\ {\mu }_{13}=\left(1-g_1\right)\left(1-g_3\right){\varphi }_S{\varphi }_I{\mu }_{av}\\ {\mu }_{33}={\left[\left(1-g_3\right){\varphi }_I\right]}^2{\mu }_{av}+{\mu }_{im}\end{array}\right\}$$

$$\left.\begin{array}{l}{C}_{12}=\left(1-c_1\right){\varphi }_S{\varphi }_F{K}_{av},C_{13}=\left(1-c_1\right)\left(1-c_3\right){\varphi }_S{\varphi }_I{K}_{av}\\ C_{23}=\left(1-c_3\right){\varphi }_I{\varphi }_F{K}_{av},K_1={\left(\left(1-c_1\right){\varphi }_S\right)}^2{K}_{av}+K_{sm}\\ K_2={\varphi }_F^2{K}_{av},K_3={\left[\left(1-c_3\right){\varphi }_I\right]}^2{K}_{av}+K_{im}\end{array}\right\}$$

$$\left.\begin{array}{l}{c}_1=K_{sm}/\left({\varphi }_S{K}_S\right)�g_1={\mu }_{sm}/\left({\varphi }_S{\mu }_S\right)\\ c_3=K_{im}/\left({\varphi }_I{K}_I\right)�g_3={\mu }_{im}/\left({\varphi }_I{\mu }_I\right)\\ K_{av}^{-1}=\left(1-c_1\right){\varphi }_F/K_S+{\varphi }_F/K_F+\left(1-c_3\right){\varphi }_I/K_I\\ {\mu }_{av}^{-1}=\left(1-g_1\right){\varphi }_S/{\mu }_S+{\varphi }_F/(2\omega {\eta }_F)+\left(1-g_3\right){\varphi }_I/{\mu }_I\end{array}\right\}$$

$$\left.\begin{array}{l}{K}_{im}={\varphi }_I{K}_I/\left[1+\varpi \left(1-{\varphi }_I\right)\right]\\ {\mu }_{im}={\varphi }_I{\mu }_I/\left[1+\varpi \gamma \left(1-{\varphi }_I\right)\right]\\ K_{sm}=\left(1-{\varphi }_F-\epsilon {\varphi }_I\right)K_S/\left[1+\varpi \left({\varphi }_F+\epsilon {\varphi }_I\right)\right]\\ {\mu }_{sm}=\left(1-{\varphi }_F-\epsilon {\varphi }_I\right){\mu }_S/\left[1+\varpi \gamma \left({\varphi }_F+\epsilon {\varphi }_I\right)\right]\\ \gamma =(1+2\varpi )/(1+\varpi ),\varpi =(1+v)/2(1-2v)\end{array}\right\}$$

$$\left.\begin{array}{l}{R}_{11}={\left[\left(1-c_1\right){\varphi }_S\right]}^2{K}_{av}+K_{sm}+4{\mu }_{11}/3\\ R_{12}=\left(1-c_1\right){\varphi }_S{\varphi }_F{K}_{av}\\ R_{22}=K_2={\varphi }_F^2{K}_{av}�R_{23}=\left(1-c_3\right){\varphi }_I{\varphi }_F{K}_{av}\\ R_{33}={\left[\left(1-c_3\right){\varphi }_I\right]}^2{K}_{av}+K_{im}+4{\mu }_{33}/3\\ R_{13}=\left(1-c_1\right)\left(1-c_3\right){\varphi }_S{\varphi }_I{K}_{av}+2{\mu }_{13}/3\end{array}\right\}$$

where $$K_S$$, $$K_F$$, $$K_I$$ are three-phase bulk modulus; $$K_{sm}$$, $$K_{im}$$ and $${\mu }_{sm}$$, $${\mu }_{im}$$ are the bulk modulus and shear modulus of soil particle skeleton and ice particle skeleton, respectively. $$K_{av}$$ and $${\mu }_{av}$$ are the undrained bulk modulus and shear modulus of the three-phase medium; $$\epsilon$$ is the contact parameter; $$c_1$$, $$c_3$$ and $$g_1$$, $$g_3$$ are the consolidation coefficients of soil skeleton and ice skeleton, respectively. $${\mu }_S$$ and $${\mu }_I$$ are the shear modulus of soil particles and ice particles, respectively. $$\varpi$$ is the cementation parameter; $$v$$ is Poisson's ratio, $${\varphi }_a(a=S,F,I)$$ denote the volume fraction of phase a.

Coefficients ${\rho }_{ij}$ in Eq. (3)

$$\left.\begin{array}{l}{\rho }_{11}=a_{13}{\varphi }_S{\rho }_S+\left(a_{12}-1\right){\varphi }_F{\rho }_F+\left(a_{31}-1\right){\varphi }_I{\rho }_I\\ {\rho }_{22}=\left(a_{12}+a_{13}-1\right){\varphi }_F{\rho }_F\\ {\rho }_{33}=\left(a_{13}-1\right){\varphi }_S{\rho }_S+\left(a_{23}-1\right){\varphi }_F{\rho }_F+a_{31}{\varphi }_I{\rho }_I\\ {\rho }_{12}=-\left(a_{12}-1\right){\varphi }_F{\rho }_F,{\rho }_{23}=-\left(a_{23}-1\right){\varphi }_F{\rho }_F\\ {\rho }_{13}=-\left(a_{13}-1\right){\varphi }_S{\rho }_S-\left(a_{31}-1\right){\varphi }_I{\rho }_I\end{array}\right\}$$

$$\left.\begin{array}{l}{a}_{12}=r_{12}\frac{{\varphi }_S\left({\varphi }_F{\rho }_F+{\varphi }_I{\rho }_I\right)}{{\varphi }_F{\rho }_F\left({\varphi }_F+{\varphi }_I\right)}+\mathrm{1,}{a}_{23}=r_{23}\frac{{\varphi }_I\left({\varphi }_F{\rho }_F+{\varphi }_S{\rho }_S\right)}{{\varphi }_F{\rho }_F\left({\varphi }_F+{\varphi }_S\right)}+1\\ a_{13}=r_{13}\frac{{\varphi }_I\left({\varphi }_S{\rho }_S+{\varphi }_I{\rho }_I\right)}{{\varphi }_S{\rho }_S\left({\varphi }_S+{\varphi }_I\right)}+\mathrm{1,}{a}_{31}=r_{31}\frac{{\varphi }_S\left({\varphi }_S{\rho }_S+{\varphi }_I{\rho }_I\right)}{{\varphi }_I{\rho }_I\left({\varphi }_S+{\varphi }_I\right)}+1\end{array}\right\}$$

where $a_{ij}$ are the tortuosity, which represents the tortuosity of j relative to i phase; $r_{ij}$ characterize the microscopic characteristics of pores, for spherical particles $r_{ij}$=0.5; ${\rho }_a(a=S,F,I)$ denote the density of phase a respectively.

Coefficients $$b_{12}$$, $$b_{23}$$, $$b_{13}$$ in Eq.(3)

$$\left.\begin{array}{l}{b}_{12}={\eta }_F{\varphi }_F^2/{\kappa }_S�b_{23}={\eta }_F{\varphi }_F^2/{\kappa }_I\\ b_{13}=b_{13}^0{\left({\varphi }_I{\varphi }_S\right)}^2�{\kappa }_S={\kappa }_{S0}{S}_r^3\\ {\kappa }_I={\kappa }_{I0}{\varphi }^3/\left[\left(1-S_r^2\right){\left(1-\varphi \right)}^3\right]\end{array}\right\}$$

where ${\eta }_F$ is the fluid dynamic viscosity coefficient; ${\kappa }_S$ and ${\kappa }_I$ are the dynamic permeability coefficients of saturated frozen soil particle skeleton and ice skeleton, respectively; ${\kappa }_{S0}$ and ${\kappa }_{I0}$ are the reference values of dynamic permeability coefficient and ice permeability coefficient of saturated soil, respectively. $b_{13}^0$ is the reference value of the viscosity coefficient. All parameters in ^{Appendix 1} are described in detail in Reference (Qiu et al. 2018).

Appendix 2

Coefficients $${\beta }_1$$, $${\beta }_2$$, $${\beta }_3$$, $${\beta }_4$$ in Eq.(21)

$${\beta }_1=-R_{11}{R}_{22}{R}_{33}+R_{11}{R}_{23}^2+R_{12}^2{R}_{33}-2R_{12}{R}_{13}{R}_{23}+R_{13}^2{R}_{22}$$,

$$\begin{array}{l}{\beta }_2=-A_{11}{R}_{22}{R}_{33}+A_{11}{R}_{23}^2+2A_{12}{R}_{12}{R}_{33}-2A_{12}{R}_{13}{R}_{23}-2A_{13}{R}_{12}{R}_{23}+2A_{13}{R}_{13}{R}_{22}-A_{22}{R}_{11}{R}_{33}+A_{22}{R}_{13}^2+\\ 2A_{23}{R}_{11}{R}_{23}-2A_{23}{R}_{12}{R}_{13}-A_{33}{R}_{11}{R}_{22}+A_{33}{R}_{12}^2\end{array}$$,

$$\begin{array}{l}{\beta }_3=-A_{11}{A}_{22}{R}_{33}+2A_{11}{A}_{23}{R}_{22}-A_{11}{A}_{33}{R}_{22}+A_{12}^2{R}_{33}-2A_{12}{A}_{13}{R}_{23}-2A_{12}{A}_{23}{R}_{13}+2A_{12}{A}_{33}{R}_{12}+A_{13}^2{R}_{22}+\\ 2A_{13}{A}_{22}{R}_{13}-2A_{13}{A}_{23}{R}_{12}-A_{22}{A}_{33}{R}_{11}+A_{23}^2{R}_{11}\end{array}$$,

$${\beta }_4=-A_{11}{A}_{22}{A}_{33}+A_{11}{A}_{23}^2+A_{12}^2{A}_{33}-2A_{12}{A}_{13}{A}_{23}+A_{13}^2{A}_{22}$$.

Coefficients $${\beta }_5$$, $${\beta }_6$$, $${\beta }_7$$ in Eq.(27)

$${\beta }_5=-{\mu }_{11}{\mu }_{33}{B}_{22}+{\mu }_{13}^2{B}_{22}$$,

$${\beta }_6=-{\mu }_{11}{B}_{22}{B}_{33}+{\mu }_{11}{B}_{23}{B}_{32}-{\mu }_{13}{B}_{12}{B}_{23}+{\mu }_{13}{B}_{13}{B}_{22}-{\mu }_{13}{B}_{21}{B}_{32}+{\mu }_{13}{B}_{22}{B}_{31}-{\mu }_{33}{B}_{11}{B}_{22}+{\mu }_{33}{B}_{12}{B}_{21}$$,

$${\beta }_7=-B_{11}{B}_{22}{B}_{33}+B_{11}{B}_{23}{B}_{32}+B_{12}{B}_{21}{B}_{33}-B_{12}{B}_{23}{B}_{31}-B_{13}{B}_{21}{B}_{32}+B_{13}{B}_{22}{B}_{31}$$.

Coefficients ${\delta }_{pn}^F$, ${\delta }_{pn}^I\left(n=\mathrm{1,}\mathrm{2,}3\right)$ in Eq. (23)

$${\delta }_{pn}^F=\frac{\left(R_{11}{R}_{23}-R_{12}{R}_{13}\right)d_n^2+\left(A_{11}{R}_{23}-A_{12}{R}_{13}-A_{13}{R}_{12}+A_{23}{R}_{11}\right)d_n+A_{11}{A}_{23}-A_{12}{A}_{13}}{\left(-R_{12}{R}_{23}+R_{13}{R}_{22}\right)d_n^2+\left(-A_{12}{R}_{23}+A_{13}{R}_{22}+A_{22}{R}_{13}-A_{23}{R}_{12}\right)d_n-A_{12}{A}_{23}+A_{13}{A}_{22}}$$,

$${\delta }_{pn}^I=\frac{\left(-R_{11}{R}_{22}+R_{12}^2\right)d_n^2+\left(-A_{11}{R}_{22}+2A_{12}{R}_{12}-A_{22}{R}_{11}\right)d_n-A_{11}{A}_{22}+A_{12}^2}{\left(-R_{12}{R}_{23}+R_{13}{R}_{22}\right)d_n^2+\left(-A_{12}{R}_{23}+A_{13}{R}_{22}+A_{22}{R}_{13}-A_{23}{R}_{12}\right)d_n-A_{12}{A}_{23}+A_{13}{A}_{22}}$$.

Coefficients $${\delta }_{sn}^F$$, $${\delta }_{sn}^F\left(n=\mathrm{1,}2\right)$$ in Eq.(29)

$${\delta }_{sn}^F=\frac{\left({\mu }_{11}{t}_n+B_{11}\right)B_{23}-B_{21}\left({\mu }_{13}{t}_n+B_{13}\right)}{-B_{23}{B}_{12}+B_{22}\left({\mu }_{13}{t}_n+B_{13}\right)}$$,$${\delta }_{sn}^I=\frac{\left(-{\mu }_{11}{t}_n+B_{11}\right)B_{22}-B_{12}{B}_{21}}{-B_{23}{B}_{12}+B_{22}\left({\mu }_{13}{t}_n+B_{13}\right)}$$.

Appendix 3

The elements of matrix [S] in Eq. (35) are:

$S_{0101}=S_{0106}=C_1$, $S_{0102}=S_{0107}=C_2$, $S_{0103}=S_{0108}=C_3$, $S_{0104}=-S_{0109}=O_1$, $S_{0105}=-S_{0110}=O_2$, $S_{0201}=-S_{0206}=F_1$, $S_{0202}=-S_{0207}=F_2$, $S_{0203}=-S_{0208}=F_3$, $S_{0204}=S_{0209}=G_1$, $S_{0205}=S_{0210}=G_2$, $S_{0301}=S_{0306}=H_1$, $S_{0302}=S_{0307}=H_2$, $S_{0303}=S_{0308}=H_3$, $S_{0304}=S_{0305}=S_{0309}=S_{0310}=0$, $S_{0401}=S_{0406}=f_{I1}$, $S_{0402}=S_{0407}=f_{I2}$, $S_{0403}=S_{0408}=f_{I3}$, $S_{0404}=-S_{0409}=I_1$, $S_{0405}=-S_{0410}=I_2$, $S_{0501}=-S_{0506}=J_1$, $S_{0502}=-S_{0507}=J_2$, $S_{0503}=-S_{0508}=J_3$, $S_{0504}=S_{0509}=P_1$, $S_{0505}=S_{0510}=P_2$, $$S_{0601}=S_{0606}=i\xi$$, $$S_{0602}=S_{0607}=i\xi$$, $$S_{0603}=S_{0608}=i\xi$$, $$S_{0604}=-S_{0609}=r_1$$, $$S_{0605}=-S_{0610}=r_2$$, $S_{0701}=-S_{0706}=-{\lambda }_1$, $S_{0702}=-S_{0707}=-{\lambda }_2$, $S_{0703}=-S_{0708}=-{\lambda }_3$, $S_{0704}=S_{0705}=S_{0709}=S_{0710}=i\xi$, $S_{0801}=-S_{0806}=-{\lambda }_1{\delta }_{p1}^F$, $S_{0802}=-S_{0807}=-{\lambda }_2{\delta }_{p2}^F$, $S_{0803}=-S_{0808}=-{\lambda }_3{\delta }_{p3}^F$, $S_{0804}=S_{0809}=i\xi {\delta }_{s1}^F$, $S_{0805}=S_{0810}=i\xi {\delta }_{s2}^F$, $S_{0901}=S_{0906}=i\xi {\delta }_{p1}^I$, $S_{0902}=S_{0907}=i\xi {\delta }_{p2}^I$, $S_{0903}=S_{0908}=i\xi {\delta }_{p3}^I$, $S_{0904}=-S_{0909}=r_1{\delta }_{s1}^I$, $S_{0905}=-S_{0910}=r_2{\delta }_{s2}^I$, $S_{1001}=-S_{1006}=-{\lambda }_1{\delta }_{p1}^I$, $S_{1002}=-S_{1007}=-{\lambda }_2{\delta }_{p2}^I$, $S_{1003}=-S_{1008}=-{\lambda }_3{\delta }_{p3}^I$, $S_{1004}=S_{1009}=i\xi {\delta }_{s1}^I$, $S_{1005}=S_{1010}=i\xi {\delta }_{s2}^I$.

where $C_n=-\text{i}\xi {\lambda }_n\left(2{\mu }_{11}+{\mu }_{13}{\delta }_{pn}^I\right)\left(n=\mathrm{1,2,3}\right)$ $O_n=-\left(r_n^2+{\xi }^2\right)\left({\mu }_{11}+\frac{1}{2}{\mu }_{13}{\delta }_{sn}^I\right)\left(n=\mathrm{1,2}\right)$, $$F_n=-i\xi {\lambda }_n\left(2{\mu }_{11}+{\mu }_{13}{\delta }_{pn}^I\right)\left(n=\mathrm{1,}\mathrm{2,}3\right)$$, $G_n=-\left(r_n{}^2+{\xi }^2\right)\left({\mu }_{11}+\frac{1}{2}{\mu }_{13}{\delta }_{sn}^I\right)\left(n=\mathrm{1,}2\right)$, $H_n=\left(R_{12}+R_{22}{\delta }_{pn}^F+R_{23}{\delta }_{pn}^I\right)\left({\lambda }_n^2-{\xi }^2\right)\left(n=\mathrm{1,}\mathrm{2,}3\right)$, $I_n=-i\xi r_n\left({\mu }_{13}+2{\mu }_{33}{\delta }_{sn}^I\right)\left(n=\mathrm{1,}2\right)$, $J_n=-i\xi {\lambda }_n\left(2{\mu }_{33}{\delta }_{pn}^I+{\mu }_{13}\right)\left(n=\mathrm{1,}\mathrm{2,}3\right)$, $P_n=-\left(r_n^2+{\xi }^2\right)\left({\mu }_{33}{\delta }_{sn}^I+\frac{1}{2}{\mu }_{13}\right)\left(n=\mathrm{1,}2\right)$.

Appendix 4

The elements in the support stiffness matrix $[M^1]$ are:

$M_{0101}^1=M_{0106}^1=C_1$, $M_{0102}^1=M_{0107}^1=C_2$, $M_{0103}^1=M_{0108}^1=C_3$, $M_{0104}^1=-M_{0109}^1=O_1^1$,$M_{0105}^1=-M_{0110}^1=O_2^1$, $M_{0201}^1=-M_{0206}^1=F_1^1$, $M_{0202}^1=-M_{0207}^1=F_2^1$, $M_{0203}^1=-M_{0208}^1=F_3^1$,$M_{0204}^1=M_{0209}^1=G_1^1$, $M_{0205}^1=M_{0210}^1=G_2^1$, $M_{0301}^1=M_{0306}^1=H_1^1$, $M_{0302}^1=M_{0307}^1=H_2^1$, $M_{0303}^1=M_{0308}^1=H_3^1$,$M_{0304}^1=M_{0305}^1=M_{0309}^1=M_{0310}^1=0$, $M_{0401}^1=M_{0406}^1=f_{I1}^1$, $M_{0402}^1=M_{0407}^1=f_{I2}^1$, $M_{0403}^1=M_{0408}^1=f_{I3}^1$,$M_{0404}^1=-M_{0409}^1=I_1^1$, $M_{0405}^1=-M_{0410}^1=I_2^1$, $$L_{0501}^1=-L_{0506}^1=J_1^1$$, $$M_{0502}^1=-M_{0507}^1=J_2^1$$,$$M_{0503}^1=-M_{0508}^1=J_3^1$$, $$M_{0504}^1=M_{0509}^1=P_1^1$$, $$M_{0505}^1=M_{0510}^1=P_2^1$$.

The elements in the support stiffness matrix $[M^N]$ are:

$M_{0101}^N=M_{0102}^N=M_{0103}^N=M_{0106}^N=M_{0107}^N=M_{0108}^N=i\xi$,$M_{0104}^N=-M_{0109}^N=r_1^N$,$M_{0105}^N=-M_{0110}^N=r_2^N$,$M_{0201}^N=-M_{0206}^N=-{\lambda }_1^N$,$M_{0202}^N=-M_{0207}^N=-{\lambda }_2^N$,$M_{0207}^N=-M_{0208}^N=-{\lambda }_3^N$,$M_{0204}^N=M_{0205}^N=M_{0209}^N=M_{0210}^N=i\xi$,$M_{0301}^N=-M_{0306}^N=-{\lambda }_1^N{\delta }_{p1}^{F^N}$,$M_{0302}^N=-M_{0307}^N=-{\lambda }_2^N{\delta }_{p2}^{F^N}$,$M_{0303}^N=-M_{0308}^N=-{\lambda }_3^N{\delta }_{p3}^{F^N}$,$M_{0304}^N=M_{0309}^N=i\xi {\delta }_{s1}^{F^N}$,$M_{0305}^N=M_{0310}^N=i\xi {\delta }_{s2}^{F^N}$,$M_{0401}^N=M_{0406}^N=i\xi {\delta }_{p1}^{I^N}$,$M_{0402}^N=M_{0407}^N=i\xi {\delta }_{p2}^{I^N}$,$M_{0403}^N=M_{0408}^N=i\xi {\delta }_{p3}^{I^N}$,$M_{0404}^N=M_{0409}^N=r_1^N{\delta }_{s1}^{I^N}$,$M_{0405}^N=M_{0410}^N=r_2^N{\delta }_{s1}^{I^N}$,$M_{0501}^N=-M_{0506}^N=-{\lambda }_1^N{\delta }_{p1}^{I^N}$,$M_{0502}^N=-M_{0507}^N=-{\lambda }_2^N{\delta }_{p2}^{I^N}$,$M_{0503}^N=-M_{0508}^N=-{\lambda }_3^N{\delta }_{p3}^{I^N}$,$M_{0504}^N=M_{0509}^N=i\xi {\delta }_{s1}^{I^N}$,$M_{0505}^N=M_{0510}^N=i\xi {\delta }_{s2}^{I^N}$.

Acknowledgments

The authors gratefully acknowledge the financial support of the Chinese Natural Science Foundation (Grant No. 52168053), the Qinghai Province Science and Technology Department Project (No.2021-ZJ-943Q), and the Project funded by China Postdoctoral Science Foundation (No.2022MD723837). The authors are also appreciating the helpful suggestions and comments of the reviewers.

Reference

• Ai, Z. Y., Ren, G. P. (2016) Dynamic analysis of a transversely isotropic multilayered half-plane subjected to a moving load. Soil Dyn Earthq Eng 83, 162-166.
• Ai, Z.Y., Gu, G.L., Zhao, Y.Z., et al. (2021). An exact solution to layered transversely isotropic poroelastic media under vertical rectangular moving loads. Comput Geotech 138, 104314.
• Ba, Z., Kang, Z., Lee, W.V. (2017) Plane strain dynamic responses of a multi-layered transversely isotropic saturated half-space. Int J Eng Sci 11955-77.
• Beskou, N.D., Theodorakopoulos, D.D. (2011). Dynamic effects of moving loads on road pavements: a review. Soil Dyn Earthq Eng 31(4), 547-567.
• Biot, M.A. (1956a). Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range. J Acoust Soc Am 28(2), 168-178.
• Biot, M.A. (1956b). Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range. J Acoust Soc Am 28(2), 179-191.
• Biot, M.A. (1962). Mechanics of deformation and acoustic propagation in porous media. J Appl Phys 33(4):1482 - 1498.
• Cao, L., Zhou, B., Li, Q., et al. (2019). Vertically dynamic response of an end-bearing pile embedded in a frozen saturated porous medium under impact loading. Shock Vib 2019, 1-18.
• Carcione, M.J., Gurevich, B., Cavallini, F. (2000). A generalized Biot-Gassmann model for the acoustic properties of shaley sandstones. Geophys Prospect 48(3), 539-557.
• Carcione, M.J., Santos, E. J., Ravazzoli, L. C., et al. (2003). Wave simulation in partially frozen porous media with fractal freezing conditions. J Appl Phys 94(12), 7839-7847.
• Carcione, M.J., Seriani, G. (2001). Wave simulation in frozen porous media. J Comput Phys 170(2), 676-695.
• Chen, C., Wang, Z.Q., Wu, W.B., et al. (2023) Semi-analytical solution for the vertical vibration of a single pile embedded in a frozen poroelastic half-space. Appl Sci 13(3), 1508-1523.
• Dong Z.J., Quan W.W., Ma X.Y., et al. (2021) Wave Propagation Approach for Elastic Transient Responses of Transversely Isotropic Asphalt Pavement under an Impact Load: A Semi analytical Solution. J Transp Eng B-Pave 147(3), 04021021.
• Eason, G. (1965). The stresses produced in a semi-infinite solid by a moving surface force. Int J Eng Sci 2(6), 581-609.
• Fan, H.S., Zhang, J.H., Zheng, J.L. (2022).Dynamic response of a multi-layered pavement structure with subgrade modulus varying with depth subjected to a moving load. Soil Dyn Earthq Eng 160, 107358.
• Jamin, P., OhSung, K., Luigi, S.D. (2021). Influence of seasonal soil temperature variation and global warming on the seismic response of frozen soils in permafrost regions. Earthq Eng Struct Dyn 50(14), 1-17.
• Jiang, H.P., Ma, Q, Shao, S.G., et al. (2023a). Characteristic of energy transmission of plane-S-wave at interface between elastic medium and saturated frozen soil medium. Chin J Rock Mech Eng 42(04):976-992.
• Jiang, H.P., Ma, Q., Cao, Y.P. (2023b). Study on the reflection and transmission of P wave on the interface between elastic medium and saturated frozen soil medium. Rock Soil Mech 44(03), 916-929.
• Leclaire, P.H. (1994) Extension of Biot's theory of wave propagation to frozen porous media. J Acoust Soc Am 96(6), 3753-3768.
• Leclaire, P.H., Frédéric, Cohen-Tenoudji., Aguirre-Puente, J. (1998). Observation of two longitudinal and two transverse waves in a frozen porous medium. J Acoust Soc Am 97(4), 2052-2055.
• Li, H.Y., Yang, S.P., Liu, J., et al. (2015). Dynamic response in multilayered viscoelastic medium generated by moving distributed loads. Eng Mech 32(01), 120-127.
• Li, Q., Shu, W., Cao, L., et al. (2019). Vertical vibration of a single pile embedded in a frozen saturated soil layer. Soil Dyn Earthq Eng 122, 185-195.
• Li, Y.C., Feng, S.J., Chen, H.X., et al. (2020). Dynamic response of a stratified transversely isotropic half-space with a poroelastic interlayer due to a buried moving source. Appl Math Model 82.
• Liu, K.F., Zhang, Z.Q., Pan, E.N. (2022). Dynamic response of a transversely isotropic and multilayered poroelastic medium subjected to a moving load. Soil Dyn Earthq Eng 155, 107154.
• Liu, S.W., Zhang, J.M. (2012) Review on physic-mechanical properties of warm frozen soil. Journal of glaciology and geocryology 34(1), 120-129.
• Ma, Q., Huang, Y.Y., Zhou, F.X., et al. (2023a). Dynamic response study of layered unsaturated foundation under moving load based on TRM method. Eng Mech.
• Ma, Q., Jiang, H.P., Zhou F.X. (2023b) Reflection and transmission of plane harmonic P wave at planar interface between elastic medium and frozen poroelastic medium, Geophys J Int 234(2), 948-971.
• Qiu, H.M. (2019). Research on propagation characteristics of elastic waves in Biot-type three-phase medium. Zhejiang: Zhejiang University, China.
• Qiu, H.M., Xia, T.D., Zheng, Q.Q., et al. (2018). Parametric studies of body waves propagation in saturated frozen soil. Rock Soil Mech 39(11), 4053-4062.
• Sheng, X., Jones, C.J.C., Thompson, D.J. (2004). A theoretical study on the influence of the track on train-induced ground vibration. J Sound Vib 272(3-5), 909-936.
• Steeb, H., Kurzeja, P.S., Schmalholz, S.M. (2014). Wave propagation in unsaturated porous media. Acta Mech 225(8), 2435-2448.
• Tabatabaie, Y.J., Valliappan, S., Zhao, C.B. (1994). Analytical and numerical solutions for wave propagation in water-saturated porous layered half-space. Soil Dyn Earthq Eng 13(4), 249-257.
• Tang, C.X., Lu, Z., Duan, Y.H., Yao, H.L. (2020). Dynamic responses of the pavement-unsaturated poroelastic ground system to a moving traffic load. Transp. Geotech 25, 100404.
• Technical Standard of Highway Engineering: JTG B01—2014[S]. Beijing: China Communications Press, 2015.
• Theodorakopoulos, D.D., Chassiakos, A.P., Beskos, D.E. (2004). Dynamic effects of moving load on a poroelastic soil medium by an approximate method. Int J Solid Struct 41, 1801-1822.
• Thomson, W.T. (1950). Transmission of elastic waves through a stratified soil medium. Appl Phys 21, 89-93.
• Wang, Y., Gao, G.Y., Yang, J., et al. (2015). The influence of the degree of saturation on dynamic response of a cylindrical lined cavity in a nearly saturated medium. Soil Dyn Earthq Eng 71(8), 27-30.
• Yang, S., Jia, M.C. (2023) Analytical layer-element solution for layered transversely isotropic saturated media subjected to rectangular moving loads. Soil Dyn Earthq Eng 171, 107877-101890.
• Yang, Y.B., Hung, H.H., Chang, D.W. (2003). Train-induced wave propagation in layered soils using finite/infinite element simulation. Soil Dyn. Earthq Eng 23 (4), 263-278.
• Ye, Z., Ai, Z.Y. (2021). Poroelastodynamic response of layered unsaturated media in the vicinity of a moving harmonic load. Comput Geotech 138:104358
• Zhang, M., Wang, X.H., Yang, G.C., et al. (2014a). Solution of dynamic Green's function for unsaturated soil under internal excitation. Soil Dyn Earthq Eng 64(64), 63-84.
• Zhang, Y., Zhang, S.X., Xia, J.H. (2014b). Transient responses of porous media under moving surface impulses. Int J Solid Struct 51, 660-72.
• Zhou, B. (2020). Study on ground motion and Rayleigh wave propagation characteristics in frozen soil area. Master's Thesis, Zhoushan: Zhejiang Ocean University, China.
• Zhou, F.X., Lai, Y.M. (2011). Propagation characteristics of elastic waves in saturated frozen soil. Rock Soil Mech 32(09), 2669-2674.

Edited by

Editor: Rogério José Marczak

Publication Dates

• Publication in this collection
  12 Feb 2024
• Date of issue
  2024

History

• Received
  07 Dec 2023
• Reviewed
  18 Jan 2024
• Accepted
  24 Jan 2024