Following Biot (1962) and Zienkiewicz et al. (1980), the equations of motion for the poroelastic medium can be expressed in terms of Cartesian coordinate system as

where the subscript index following a comma and dot above a symbol indicate the derivative with respect to a spatial coordinate and time, respectively; σ_ij is the total stress tensor; p is the fluid pressure; ϕ is the porosity; w_i =ϕ(U_i -u_i ) is the average displacement component of the fluid relative to the solid in which u_i and U_i are the displacement components of the solid skeleton and pore fluid, respectively; ρ=(1-ϕ)ρ_s +ϕρ_f is the density of the solid-fluid mixture with ρ_f and ρ_s being the densities of fluid and solid, respectively; k_i is the permeability coefficient (or called hydraulic conductivity coefficient) in the ith direction; g indicates the acceleration of gravity.

Following generalized Hooke’s law, the stress-strain relationship under plane strain condition can be written as

in which

where α_i (i=1, 3) and M are the Biot’s effective stress coefficients and Biot’s modulus, respectively; C ₁₁, C ₁₂, C ₁₃, C ₃₃, C ₄₄ are the elastic constants, and the relation between the elastic constants and engineering parameters is listed in Appendix A ; K_s and K_f are bulk moduli of solid skeleton and the pore fluid, respectively.

When seabed is in the unsaturated state with very small amount of gas, the following relation holds

where S_r is the degree of saturation; K_w is the true bulk modulus of elasticity of water which is generally selected as 2×10⁹Pa; P_{w 0}=ρ_fgd is the absolute pore-fluid pressure with d being the depth of the sea water.

As already mentioned, we have assumed that any layer interface of the multilayered seabed is perfectly connected. Hence, the continuity conditions on z=z_f can be written as (i=x, z; j=x, z)

Since the wave-induced response of the seabed is periodic, we express the field quantities in the form of complex variables as

where m=1, 2 represents the response caused by the first-order and the second-order waves, respectively. U ¯ x ( m ) , U ¯ z ( m ) , W ¯ x ( m ) , W ¯ z ( m ) and P ¯ ( m ) are the magnitudes of the dynamic response induced by the wave loading. To facilitate the derivation for the solution of the layered system in the following section, the stress components (e.g., σ_zz (x, z, t), τ_xz (x, z, t)) should also be expressed as the similar complex variable form. It is noted that the solution in the complex variables domain can be solved first, then taking the summation of the real part of each order results in the final solution.

Substituting Eq. (14) into Eq. (2), the following relations can be derived

where δ 1 ( m ) = ( ρ f m 2 ω 2 / ϕ + i m ω ρ f g / k x ) − 1 ;

Combining Eqs. (1), (3), (14)-(16) and eliminating w_x and w_z , we have

where D is the differential operator (i.e., D=∂/∂z); a i ( m ) (i=1-6) are the coefficients with being defined as

The potential function Φ ¯ ( m ) ( z ) is introduced, which has been successfully applied to the dynamic response of the seabed under linear wave (Li et al., 2020). The aforementioned field quantities can be expressed in terms of potential function as

It is noted that the introduced potential function automatically satisfies the first and third equations in Eq. (18). Substituting Eq. (20) into Eq. (18), the final form of the governing equation can be written as

where r i ( m ) (i=1-4) are coefficients, whose detailed expressions are given in Appendix B .

Through some algebraic calculations, the solution of Eq. (21) can be given as

where ± λ i ( m ) (i=1-3) are given in Appendix C . The values of A i ( m ) and B i ( m ) (i=1-3) are to be determined by the boundary conditions.

Substituting Eq. (22) into Eq. (20) and performing some algebraic calculations, the general solutions of displacements, pore pressure and stresses in any homogeneous layer under the action of second-order Stokes wave can be expressed as

where χ i ( m ) , φ i ( m ) , ξ i ( m ) and η i ( m ) (i=1-3) are the coefficients given in Appendix D . It is noted that the effective stresses in Eq. (24b) are derived via the relations between total stress and pore pressure.

We adopt the dual variable and position (DVP) method (Pan, 2019; Liu et al., 2022) to expand the single-layer solution to the multilayered solution. It is stated again that the layered solution corresponding to different m will be solved separately and then be superposed together to gain the complete solution. To facilitate the derivation, e^{i m ( kx - ωt )} is suppressed, and the vectors and diagonal matrices are defined as

U ( m ) ( z ) = [ U ¯ x ( m ) ( z ) , U ¯ z ( m ) ( z ) , W ¯ z ( m ) ( z ) ] T ;

E 1 ( m ) ( z ) = d i a g [ e λ 1 ( m ) z , e λ 2 ( m ) z , e λ 3 ( m ) z ] ;

where τ ¯ x z ( m ) ( z ) and σ ¯ z z ( m ) ( z ) are the stress magnitudes.

Then the single-layer solution for jth layer can be rewritten in terms of matrix form as

where [ M i j ( m ) ] is the 3×3 submatrix of 6×6 matrix M^(m) , the elements of which are given in Appendix E .

Substituting z=z_{j −1} and z=z_j into Eq. (27), the solutions for the top and bottom interfaces of jth layer can be written as

By eliminating the unknown vectors K + ( m ) , K − ( m ) and making use of DVP method, the relation for the layer j can be expressed as

where

After gaining the layer-matrix relation for layer j+1, utilizing the continuity conditions Eq. (7) and making some algebraic operations leads to the recursive relationship from jth layer to (j+1)th layer

in which

where I denotes the identity matrix.

After obtaining the recurrence relation Eq. (32), the dynamic response of the layered seabed can be solved according to the boundary conditions. For the seabed of finite thickness, the boundary conditions in Eqs. (12) and (13) can be rewritten as

where the subscript 1 in α ₃₁ denotes the layer number; P ⁽¹⁾=P ₁ and P ⁽²⁾=P ₂.

In order to solve the dynamic response at arbitrary depth z=z_a (say in layer j, see Figure 1 ), we further divide the homogeneous layer j into two sublayers, i.e., sublayers a1 and a2 on the top and bottom parts, respectively. We propagate the recursive relation Eq. (32) from layer 1 to sublayer a1 and from sublayer a2 to layer n, which yields the global matrix as

Till here, the solutions for the filed quantities corresponding to different m are gained. Then taking the summation of solutions corresponding to m=1 and 2 yields the final solutions. It should be pointed out that for the homogeneous TI poroelastic seabed, the analytical solution can be easily derived by combining the boundary conditions in Eqs. (12) and (13), and general solution in Eqs. (23) and (24). In addition, if the seabed is actually a layered half-space, we only need set the nth layer with a very large thickness to gain the half-space results.

In order to validate the reliability of the present solution, two cases corresponding to single-layer and multilayered seabed under linear wave are considered first. In the comparison, we take the amplitude of the pore pressure |p| and vertical effective stress |σ ^′ _zz | as the analyzed physical quantities. That is to say, |p| and |σ ^′ _zz | are the maximum value of p and σ ^′ _zz , respectively. The variable used for normalization is defined as p ₀=ρ_fgH/(2cosh(kd)). For the single-layer seabed (e.g., graveled seabed), the solution of Jeng and Lee (2001) and the reduced one of the present study are compared with using the following parameters: T=15 s, d=20 m, H=2 m, h=20 m, v=1/3, ϕ=0.35, S_r =0.95-1.0, ρ_s =2650 kg/m³, ρ_f =1000 kg/m³, k_x =k_z =1×10^-1 m/s, G_h =G_v =5×10⁷ N/m². Moreover, the following poroelastic properties α ₁=α ₃ = 1 and M=K_f /ϕ are used in the present solution. Figure 2 shows the comparison of the present reduced solution with that by Jeng and Lee (2001) for a homogeneous isotropic and poroelastic seabed subjected to linear wave. In order to gain the reduced solution for linear wave, we need to fix P ₂ = 0 and thus m=1. Notice that the complete Biot’s poroelastodynamic theory are considered in both studies. It is observed from Figure 2 that the solution of the present study is in good agreement with that of Jeng and Lee (2001).

We further consider a two-layer seabed under linear wave and compare the reduced results from present study with those by Hsu et al. (1995) as shown in Figure 3 . The parameters used for verification are: T=10 s, d=20 m, H=6 m, L=121.12 m, k_{z 2} = 10^-3 m/s, h ₁ = 10 m and h ₂ = 40 m. The remaining seabed parameters, except for permeability coefficient, are the same for both two layers and taken as v=1/3, ϕ=0.3, S_r =0.975, ρ_f =1000 kg/m³, ρ_s =2000 kg/m³, G_v =G_h =1×10⁷ N/m², α ₁=α ₃ = 1 and M=K_f /ϕ. It can be seen from Figure 3 that the results from the reduced solution of present study agree well with those of Hsu et al. (1995). It can be concluded from above two cases that the present solutions are applicable for both the single-layer and multilayered cases.

The present solution is further compared with the existing one for a homogeneous, isotropic and poroelastic seabed under nonlinear wave. In the comparison, the parameters are: T=10 s, d=10 m, L=240 m, H=0.08tanh(kd) m, h=50 m, v=0.333, ϕ=0.3, ρ_s =2650 kg/m³, ρ_f =1000 kg/m³, k_x =k_z =1×10^-3 m/s, G_h =G_v =1×10⁷ N/m², α ₁=α ₃ = 1 and M=K_f /ϕ. Figure 4 shows the comparison of the present reduced solution with that by Zhou et al. (2011) for a homogeneous, isotropic and poroelastic seabed subjected to second-order Stokes wave. |u_z | in Figure 4 denotes the amplitude of u_z . It can be observed from Figure 4 that the results from the present study have good agreement with those from Zhou et al. (2011).

It has been mentioned in Li et al. (2020) that the linear wave-induced dynamic response of a finite-thickness seabed is dependent on the properties of the soil, and anisotropic stiffness and permeability have significant effect on the seabed response. However, the dynamic response of the anisotropic seabed under non-linear wave is still unknown. Therefore, in the subsequent study, we will analyze the dynamic response and liquefaction behavior of single-layer and multilayered anisotropic seabed with different anisotropic parameters and degree of saturation and subjected to non-linear wave. The basic poroelastic properties and wave conditions used in the following analysis are listed in Table 1 , and it should be pointed out that all parameters are taken from this table if there is no further statement.

Figure 5 shows the comparison of the induced pressures induced by non-linear wave (i.e., second-order Stokes wave) and linear wave at the seabed surface. It can be observed that the wave loading by non-linear wave is different with that by linear wave. Compared to linear wave, the wave loading by non-linear wave shows higher/sharper wave crest and lower/flatter wave trough, and takes on more evident characteristics of asymmetric distribution. As a result, fully understanding the difference in the induced filed quantities (e.g., pore pressure, vertical effective stress) by linear and non-linear waves for different soil properties is of great importance. In the following parts, based on the input data given in Table 1 , the influence of anisotropic stiffness, anisotropic permeability, degree of saturation and stratification on the vertical distribution of maximum pore pressure |p|, vertical effective stress |σ ^′ _zz | and shear stress |τ _xz | are analyzed for both linear and non-linear waves. Unless otherwise stated, the solid and dash line denote, respectively, the response of non-linear and line waves in Figures 6 – 14 .