2D Numerical Study of the Stability of Trench under Wave Action in the Immersing Process of Tunnel Element

Abstract

The evaluation of the trench stability under the action of ocean waves is an important issue in the construction of an immersed tunnel. In this study, a two-dimensional coupling model of a wave-seabed-immersed tunnel is proposed for the dynamic responses of a trench under wave action in the immersing process of tunnel elements. The porous seabed is characterized by Biot consolidation equations. The $k-\epsilon$ model and RANS equation are adopted to achieve the flow field simulation, and the level set method (LSM) is used to capture the free surface between the water and air. The proposed numerical model is verified using the experimental data and analytical results. Then, the transient liquefaction and shear failure in the vicinity of the trench are discussed at two different conditions, namely, after the foundation groove is excavated and after the tunnel element is placed. The pore pressure amplitude on the weather side slope is demonstrated to be significantly smaller than that on the lee side slope. Also, the distribution of the surrounding flow field and pressure field change dramatically after the tunnel element is settled, leading to the significant changes of seabed stability.

1. Introduction

With the development of ground transportation, inland river transportation and sea transportation, there is an increasing demand for channels across the river and sea. An immersed tube tunnel, as a kind of underwater tunnel construction method, has been developed since the beginning of the 20th century. Its applicability and reliability have been verified through many successful cases of immersed tube tunnels around the world such as the Oresund tunnel between Denmark and Sweden, as well as the Hong Kong-Zhuhai-Macao immersed tunnel in China. By excavating a trench of fixed geometry into the soil surface, the prefabricated pipe sections are laid successively in the trench. The trench soil is usually soft and very sensitive to the contact deformation and load. Moreover, the structure is also sensitive to the uneven settlement of the foundation as the immersed tube tunnel could be regarded as nearly infinite in length direction compared with its section size. Thus, the stability of the excavated trench soil around the immersed pipe tunnel is an important issue for the safety of the tunnel.

Extensive research has been conducted to the dynamic interaction mechanism between wave/flow, seabed and structures [1,2,3,4,5,6,7,8,9,10,11]. In general, when the wave propagates over the seabed, the dynamic wave pressure will induce pore water pressure in the seabed. When the pore pressure reaches the limit value, the effective stress in the soil disappears and is always accompanied by instability or even liquefaction. According to the available literatures, the wave-induced liquefaction can be divided into two forms according to the generation mechanisms of the pore pressure. One is due to the accumulative effect of pore pressure, caused by the volume deformation of the soil under cyclic loading [12,13]. Another is induced by transient pore water pressure, which usually happens under wave troughs, with the reduction of pore pressure amplitude and phase lag [14,15,16].

In recent years, many researchers used various analytical formulas, which assume that the seabed is flat, to study the dynamic response of the porous seabed under linear waves [17,18,19]. However, this method is difficult to consider the existence of the structure. Meanwhile, a numerical simulation is also used to study the response mechanism of the wave-sea-floor structure system. Jeng and Cheng [20] used the finite difference method to study the effect of the wave on the soil around a buried pipe. Jeng et al. [6] studied the coupled response of the wave-seabed-breakwater system by using the finite element method. Zhao et al. [21] used the combined FVM-FEM (Finite Volume Model-Finite Element Model) to study the effect of the flow field on the submerged rubble mound breakwater. Kasper et al. [22] used the FEM method in a practical engineering and studied the stability of an immersed tunnel under deep water wave impact. Liao et al. [23,24] adopted a two-step numerical scheme combined CFD (Computational Fluid Dynamics) model to study the oscillatory response around slope breakwater heads.

The immersed tunnel, however, is rarely concerned in the existing literature. However, the dynamic response of seabed soil in fluid field could provide a reference for this study. The size of an immersed tube tunnel is relatively large (e.g., the one-way three-lane tunnel designed in the Hong Kong-Zhuhai-Macao bridge has a single-hole span of 14.55 m). When the tunnel section is placed in the trench, the distribution of the surrounding flow field will be inevitably changed. The flow field around the tunnel could then greatly affect the stability of the surrounding soil. However, the pressure boundary setting method, associated with the analytical solution, cannot capture the local flow state or perform the bidirectional coupling process.

In this study, a coupling model of the wave-seabed-tunnel is proposed to study the transient response of the soil in the vicinity of a trench and the flow field distribution nearby. The wave field is simulated by solving the RANS equation and $k-\epsilon$ turbulence model. The dynamic response of the soil is solved by the Biot porous elastic theory, and the deformation and displacement of the tunnel are ignored in this study. In Section 2, the wave model and the soil model are respectively verified in detail. In Section 3, the dynamic response and failure of the trench soil are analyzed; the key parameters affecting the soil instability around the tunnel are also discussed.

2. Numerical Model

The numerical model mainly includes two parts: a wave model and a seabed model with an immersed tube tunnel. The wave model is used to simulate the flow field, and the seabed model is used to study the oscillatory response in the trench around the tube tunnel in the flow field. Figure 1 describes a sketch of the 2-D wave-seabed-tunnel interaction, where $h$ is the seabed depth, $d$ is the water depth, A1 and A2 are two asymmetric points on two slopes, B1 and B2 are 1 m below the positions of two bottom corners of the tunnel and $\theta$ is the slope ratio of the trench. The proposed numerical model mainly focuses on the comparison between the case of the free field and the case after the tunnel is settled. Thus, there is no backfill material in the trench. The detailed cross-sectional dimensions of the tunnel are depicted in the following section.

2.1. Fluid Dynamic Model

In this part, the momentum source function is used to make waves by replacing the source term with the momentum source function in the momentum conservation equation [25,26]. Moreover, the level set method (LSM) is used to capture the free surface between water and air, and this free contact surface is a time-dependent variable varying with time as the wave surface changes.

The $k-\epsilon$ model and RANS equation are adopted to achieve a flow field simulation. The standard turbulence model mainly includes two transport equations and two independent variables: turbulent kinetic energy $k$ and turbulent dissipation rate $\epsilon$ [27]. Turbulence is a complex, unsteady, irregular flow with rotation, which is caused by viscous forces. When the Reynolds number is large enough and the inertial force dominates, the flow becomes turbulent. The eddy viscosity is

$${\mu }_T=\rho C_{\mu }\frac{k^2}{\epsilon }$$

(1)

where $C_{\mu }$ is a model constant and $\rho$ is the fluid density.

The transport equation for turbulent kinetic energy k is written as

$$\rho \frac{\partial k}{\partial t}+\rho u\cdot \nabla k=\nabla \cdot \left(\left(\mu +\frac{{\mu }_T}{{\sigma }_k}\right)\nabla k\right)+P_k-\rho \epsilon$$

(2)

where $u$ is the instantaneous velocity in vector notation, $\nabla$ is the Laplace operator, $t$ is time, $\mu$ is the molecular viscosity and the production term can be expressed as

$$P_k={\mu }_T\left(\nabla u:\left(\nabla u+{\left(\nabla u\right)}^T\right)-\frac{2}{3}{\left(\nabla \cdot u\right)}^2\right)-\frac{2}{3}\rho k\nabla \cdot u$$

(3)

The transport equation for dissipation rate $\epsilon$ is written as

$$\rho \frac{\partial \epsilon }{\partial t}+\rho u\cdot \nabla \epsilon =\nabla \cdot \left(\left(\mu +\frac{{\mu }_T}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right)+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}$$

(4)

The aforementioned equations include closure coefficients $C_{\mu }$, $C_{\epsilon 1}$, $C_{\epsilon 2}$, ${\sigma }_k$ and ${\sigma }_{\epsilon }$ which are 0.09, 1.44, 1.92, 1.0 and 1.3, respectively.

The stress tensor, including viscous stress and Reynolds stress, can be expressed as

$${\tau }_{ij}=\upsilon \left[\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right]-\rho {\overline{u}}_i{\overline{u}}_j$$

(5)

where $\upsilon$ is the dynamic viscosity and $-\rho {\overline{u}}_i{\overline{u}}_j$ is the Reynolds stress term; under the eddy-viscosity assumption [27,28], the Reynolds stress term can be expressed as

$$\rho {\overline{u}}_i{\overline{u}}_j=-2{\mu }_T{S}_{ij}+\frac{2}{3}\rho k{\delta }_{ij}$$

(6)

where ${\delta }_{ij}$ is Kronecker delta and $S_{ij}$ is the strain-rate tensor.

For this 2-D problem in this model, the mass conservation equation and the momentum conservation equation can be expressed as

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=0$$

(7)

$$\frac{\partial \overline{u_i}}{\partial t}+\overline{u_j}\frac{\partial \overline{u_i}}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x_j}+\frac{1}{\rho }\frac{\partial {\tau }_{ij}}{\partial x_j}+g+S_i$$

(8)

where, $\overline{u_i}$, $\overline{u_j}$ and $\overline{p}$ are the average velocity in the x and z direction and the wave pressure in the flow field, respectively; $x_i$ and $x_j$ are the cartesian coordinates that are orthogonal to each other; $t$ is time; $\rho$ is the density of the fluid; $g$ is the acceleration of gravity; ${\tau }_{ij}$ is the shear stress tensor; and $S_i=(S_x,S_z)$ is the momentum source function. In this study, the momentum source function can be expressed as

$$S_x=-g(2\beta x)\exp (-\beta x^2)\frac{D_s}{\omega }\sin (k_{w}z-\omega t)$$

(9)

$$S_z=g\exp (-\beta x^2)\frac{k{D}_s}{\omega }\cos (k_{w}z-\omega t)$$

(10)

where $k_w$ denotes the wave number; $\omega$ is the wave angular frequency; and $\beta =80/{\delta }^2/L^2$ is the width of the wave source region, in which $L$ is the wave length and $\delta$ is a parameter that indicates the wave source region. $D_s$ denotes amplitude of source region and is expressed as

$$D_s=\frac{2A({\omega }^2-{\alpha }_1\mathrm{g}{k}_w{}^4{d}^3)\cos {\theta }_w}{\omega I_1{k}_w\left[1-\alpha {(k_{w}d)}^2\right]}$$

(11)

where $A$ is the wave amplitude, $d$ is the water depth, ${\theta }_w$ is the angle between the direction of wave propagation and the horizontal line (set as zero in the model), $I_1=\sqrt{\pi /\beta }\exp (-k^2/4\beta )$, ${\alpha }_1=\alpha +1/3$ and $\alpha$ is another wave parameter expressed as

$$\alpha =\frac{z_{\alpha }}{d}\left(\frac{z_{\alpha }}{2d}+1\right), z_{\alpha }=-0.53d$$

(12)

Detailed parameters of other related parameters can be found in Wei and Kirby (1999) [25] and Choi and Yoon (2009) [26].

The LSM method is used to capture the free surface between water and air, which can be expressed as

$$\frac{\partial \varphi }{\partial t}+u\cdot \nabla \varphi =\gamma \nabla \cdot \left({\epsilon }_t\nabla \varphi -\varphi \left(1-\varphi \right)\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)$$

(13)

where $\varphi$ is the level set function; ${\epsilon }_t$ is the parameter controlling the interface thickness, and its default value is half of the maximum mesh element size in which the interface passes through; $\gamma$ is the reinitialization parameter, and the default value is 1 m/s; and $u$ is the velocity field component, and the applied velocity field transports the level set method through convection. A detailed introduction about the LSM method could be found in Olsson et al. [29], Olivier et al. [30] or Sheu et al. [31].

2.2. Seabed Model

In the seabed model, the pore pressure changes periodically due to the wave pressure obtained from the upper wave model. Two main stages during the construction of the immersed tunnel are considered. One is the state when the foundation groove excavation is completed, and another is the unfilled state when the tunnel subsidence is completed. Due to the short time between the placement of the immersed tube tunnel and backfilling, the transient response caused is considered in this paper. The 2-D seabed model is established by adopting the Biot consolidation equation, which ignores the accelerating inertia term of flow and soil particles [21,23,24].

Assuming the seabed to be homogeneous and isotropic, the permeability coefficient in all directions are the same. The mass conservation equation of pore fluid can be expressed as

$$k_s{\nabla }^{2}p-{\gamma }_w{n}_s{\beta }_s\frac{\partial p_s}{\partial t}={\gamma }_w\frac{\partial {\epsilon }_s}{\partial t}$$

(14)

where $k_s$ is the Darcy permeability coefficient, ${\nabla }^2={(\frac{\partial }{\partial x})}^2+{(\frac{\partial }{\partial z})}^2$ is the Laplace operator, ${\gamma }_w$ is the unit weight of pore water, $n_s$ is the soil porosity and ${\epsilon }_s=\frac{\partial u_s}{\partial x}+\frac{\partial v_s}{\partial z}$ is the soil elastic volume strain. The compressibility of pore water could use ${\beta }_s$ expressed as [32]

$${\beta }_s=\frac{1}{K_f}=\frac{1}{K_w}+\frac{1-S_r}{p_0}$$

(15)

where $K_f$ is the measured volume modulus of pore water, $K_w$ is the true volume modulus of water, $p_0$ is the absolute water pressure (compared to zero), $p_0={\gamma }_{w}d$ is defined in the model, $d$ is the water depth and $S_r$ is the soil saturation.

The governing equation for the porous seabed can be expressed as

$$\frac{\partial {\sigma }_{xx}^{\prime }}{\partial x}+\frac{\partial {\tau }_{xz}}{\partial z}=\frac{\partial p_s}{\partial x}$$

(16)

$$\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\sigma }_{zz}^{\prime }}{\partial z}+\rho g=\frac{\partial p_s}{\partial z}$$

(17)

where ${\sigma }_{xx}^{\prime }$ and ${\sigma }_{zz}^{\prime }$ are the horizontal and vertical effective normal stresses; ${\tau }_{xz}$ is the shearing stress; the actual density of seabed is $\rho =(1-n){\rho }_s+n{\rho }_w$, in which ${\rho }_s$ is the density of soil particles and ${\rho }_w$ is the water density.

The aforementioned stress can be derived in the poro-elastic theory; the relationship between effective stress and deformation can be expressed as

$${\sigma }_{xx}^{\prime }=\frac{2G_s}{1-2{\nu }_s}\left[\left(1-{\nu }_s\right)\frac{\partial u_s}{\partial x}+{\nu }_s\frac{\partial v_s}{\partial z}\right]$$

(18)

$${\sigma }_{zz}^{\prime }=\frac{2G_s}{1-2{\nu }_s}\left[\left(1-{\nu }_s\right)\frac{\partial v_s}{\partial z}+{\nu }_s\frac{\partial u_s}{\partial x}\right]$$

(19)

$${\tau }_{xz}={\tau }_{zx}=G_s\left(\frac{\partial v_s}{\partial x}+\frac{\partial u_s}{\partial z}\right)$$

(20)

where $G_s$ is the soil shear modulus, which can be expressed through Young’s modulus $E_s$ and Poisson’s ratio ${\nu }_s$ as $G_s=E_s/\left[2\left(1+{\nu }_s\right)\right]$.

2.3. Tunnel Model

Herein, the immersed tunnel placed in an underwater trench is considered to be an impermeable uniform elastic medium with little deformation. Based on Hooke’s law, the governing equation for the immersed tunnel can be derived:

$$G_t\left(\frac{{\partial }^2{u}_t}{\partial x^2}+\frac{{\partial }^2{u}_t}{\partial z^2}\right)+\frac{G_t}{1-2{\nu }_t}\frac{\partial }{\partial x}\left(\frac{\partial u_t}{\partial x}+\frac{\partial v_t}{\partial z}\right)=0$$

(21)

$$G_t\left(\frac{{\partial }^2{v}_t}{\partial x^2}+\frac{{\partial }^2{v}_t}{\partial z^2}\right)+\frac{G_t}{1-2{\nu }_t}\frac{\partial }{\partial z}\left(\frac{\partial u_t}{\partial x}+\frac{\partial v_t}{\partial z}\right)+{\rho }_{t}g=0$$

(22)

where $u_t$ and $v_t$ are displacements in the x direction and z direction respectively, ${\rho }_t$ is the density of the tunnel, and $G_t$ and ${\nu }_t$ represent the shear modulus and Poisson’s ratio of the immersed tunnel, respectively.

2.4. Boundary Conditions

The wave field is simulated numerically using the wave model. Under wave loading, the wave pressure is transmitted to the porous seabed with a finite thickness and the bottom of the seabed is assumed to be rigid and impermeable. The still water surface is located at z = 0, the waves propagate in the positive x direction, the vertical z-axis starts at the still water surface and upward is positive. In the wave field, H is the wave height, L is the wave length, the distance between the still water surface and seabed surface is d, and the seabed thickness is h. Herein, the oscillatory pore water pressure and seabed deformation could be obtained through solving Equations (16)–(20). Thus, the boundary conditions should be defined clearly, including the seabed conditions, tunnel surface boundary conditions and free water surface boundary conditions.

2.4.1. Seabed Boundary Conditions

At the seabed surface, the vertical effective stress and shear stress become zero, and the pore water pressure is equal to the wave pressure in a wave field.

$${\sigma }_{zz}^{\prime }=0,{\tau }_{xz}^{\prime }=0,p_s=p_b\left(\mathrm{at}\mathrm{seabed}\mathrm{surface}\right)$$

(23)

where $p_b$ is the wave pressure at the seabed surface.

The bottom of the seabed is impermeable bedrock; there is no displacement and vertical seepage at the bottom of the seabed, which is expressed as

$$u_s=v_s=0,\frac{\partial p_s}{\partial z}=0(\mathrm{at}\mathrm{the}\mathrm{bottom}\mathrm{of}\mathrm{seabed})$$

(24)

Both sides of the seabed boundary are impermeable, and there is no horizontal displacement, which is expressed as

$$u_s=0,\frac{\partial p_s}{\partial z}=0$$

(25)

2.4.2. Tunnel Boundary Conditions

The tunnel is assumed to be impermeable and elastic material with large stiffness; the pore water pressure gradient at the tunnel’s outside surface is zero:

$$\frac{\partial p_s}{\partial n}=0$$

(26)

Meanwhile, no relative displacement is assumed to occur between the tunnel and seabed.

2.4.3. Free Water Surface Boundary Conditions

Owing to the pressure continuity condition of the two-phase flow in the flow field, the pressure at the water free surface is equal to the atmospheric pressure.

2.5. Integration of Fluid Dynamic Model and Seabed Model

In this model, FEM (Finite Element Model) codes are constructed within the COMSOL Multiphysics by designing the user’s own governing equations and boundary conditions in a certain physical model. To figure out the dynamic response of seabed around the immersed tunnel under the action of waves, the RANS equation with momentum source function source terms is solved to realize a wave simulation and the PDE (partial differential equation) module is defined to realize a finite element calculation for the seabed soil. In this model, the tunnel is considered to be an impermeable, uniform, elastic medium with great stiffness; the interaction between the soil and the tunnel structure is not discussed in detail. Moreover, the dynamic interaction between the fluid and the tunnel might affect the stability of the tunnel. This paper, however, focuses on the stability of the soil trench under wave action in the immersing process of the tunnel element. Thus, the interaction between the fluid and the tunnel is also not emphasized.

First, MUMPS (a multifrontal massively parallel sparse direct solver) is adopted to solve the governing equations for the wave field and seabed, and the pressure fields and displacement fields can be obtained. Then, PARDISO (a nested dissection multithreaded solver) is used to get the turbulence energy and turbulent dissipation rate due to the large amount of computation in solving these two turbulence variables. The whole wave field is simulated to realize the flow field and wave pressure field simulation around the seabed and immersed tunnel. At last, two sub-modules (i.e., the wave model and seabed model) are coupled to achieve physical field coupling. Therefore, in every time step, the soil response under the wave field could give feedback to the wave sub-model.

The meshes of the model are divided into two types. The meshes are divided into triangles automatically in the wave flow field with a soil and tunnel part. Meanwhile, the automatic subdivision method is used to match the grids well with physical fields. The meshes around the foundation trench of the tunnel are also refined, which could make the simulation results of the flow field more reliable. The moving mesh method is used to capture the interface between the free surface and air, and the Lagrange smooth type is used in the freely deformed mesh area. The wave surface at 1 m outside the source region is selected for the comparison with the analytical solution, as shown in Figure 2.

3. Model Validation and Numerical Results

3.1. Model Validation

In this section, the calculated results of the proposed model are compared with the experimental values or analytical solutions from previous studies. At first, the wave verification is carried out to ensure the accuracy of the simulated waves. The wave parameters in Liu et al. [33] are adopted: wave height $H = 0.143 \mathrm{m}$, wave period $T = 1.75 \mathrm{s}$, water depth $d = 0.533 \mathrm{m}$ and wave length $L = 3.53 \mathrm{m}$. The wave surface at 1 m outside the source region is selected for the comparison with the analytical solution, as shown in Figure 2. The results in the present wave model overestimate the wave surface at first; however, the overall result is basically consistent with the results. Meanwhile, the wave motion pattern at the free wave surface is compared with the analytical solution, as shown in Figure 3. The result of the present model tends to be consistent with the analytical solution. As a whole, the results of the wave simulation agree well with the expected analytical solution.

Then, the dynamic response in the seabed is compared with the analytical solutions [34] and experimental data [33]. Figure 4 shows the variation of the maximum pore pressure ($p_s/p_0$) with soil depth ($z/h$), in which $p_s$ denotes the wave-induced maximum pore pressure and $p_0$ is the pressure amplitude at the seabed surface. The soil parameters are given in

Table 1. The soil parameters for the present study.

shear modulus ($G$)	1.27 × 107 N/m2
poison’s ratio ($\upsilon$)	0.3
soil permeability ($k_f$)	1.8 × 10−4 m/s
soil porosity ($n$)	0.425
saturation degree ($S_r$)	0.995
seabed thickness ($h$)	1.8 m

. It is shown that the simulated result from the present numerical model agrees well with the analytical solution in Hsu and Jeng [34], while the gap between the present model and the experimental data becomes more significant near the bottom of the soil. However, this paper focuses on the dynamic response of the soil near the seabed surface which has significant influence on the stability of the trench. Therefore, the reliability of the proposed model in a soil surface layer could make the following analysis realized with good efficacy.

3.2. Consolidation of the Seabed

In a natural marine environment, the seabed soil will reach a new state of consolidation after the construction of the structure, with the dissipation of the pore pressure, the increase of effective stress and the settlement of the foundation. The stress distribution around marine structures may be significantly affected, so the initial consolidation state and stress distribution should be clearly defined. In this section, the initial consolidation of the porous seabed is predetermined under dead weight and tunnel load. The tunnel size and material parameters are given in Figure 5 and

Table 2. The input data of a standard case for parametric study.

Wave Parameters	Value	Unit
wave height ($H$)	2	m
wave period ($T$)	8	s
wave length ($L$)	83.4	m
water depth ($d$)	16	m
Soil Parameters
seabed thickness ($h$)	30.5	m
shear modulus ($G$)	5 × 106	N/m2
soil porosity ($n$)	0.45	-
poison’s ratio ($\upsilon$)	0.27	-
elastic modulus ($E$)	3 × 107	N/m2
soil permeability ($k_f$)	10−6	m/s
saturation degree ($S_r$)	0.975	-
density of soil grain (${\rho }_s$)	2650	kg/m3
internal cohesion ($c$)	0	kPa
internal friction angle ($\varphi$)	30	deg
Water Parameters
shear modulus ($G$)	2 × 109	N/m2
density of water (${\rho }_w$)	986	kg/m3
Tunnel Parameters
elastic modulus ($E_t$)	3.5 × 1010	N/m2
poison’s ratio (${\upsilon }_t$)	0.18	-
density of tunnel (${\rho }_t$)	2700	kg/m3

.

When the immersed tunnel is arranged, its sinking velocity is basically stable by controlling the irrigation amount; thus, the self-weight of the tunnel is basically equivalent to the buoyancy. The effect of the tunnel placement on the internal force of the seabed soil is not considered in this model. The initial effective stress distribution is shown in Figure 6. This pre-consolidation is considered as the initial circumstance in the following analysis.

3.3. Dynamic Responses of the Seabed

Figure 7 is given to compare the distributions of pore pressure and various stresses in two cases: t = 29 s after the foundation groove is excavated and t = 29 s after the tunnel element is placed. The horizontal/vertical stresses and shear stress, as well as the pore pressure, are respectively depicted. These response variables are nondimensionalized using the wave-induced pressure at the seabed surface. It is shown that the horizontal and vertical stresses are right below the wave crest and wave trough, and the effective horizontal stress under the wave crest is negative while the effective vertical stress is positive. The shear stress is shown to be concentrated at the joint where the crests and troughs meet (i.e., the position of wave node where the wave amplitude is zero). From the comparison between the two cases, the horizontal and vertical effective stresses in the lower part of the tunnel tend to be concentrated. The negative pore pressure region on the right side of the tunnel expands, which is consistent with the expansion of a liquefaction trend which will be discussed in the subsequent section.

Figure 8 shows the variation of pore pressure with time at 1 m below the middle point in both the lee (left) side and weather (right) side slopes. Due to the existence of the trench, the pore pressure amplitude on the weather side slope is significantly smaller than that on the lee side slope. By comparing Figure 8a with Figure 8b, the amplitude of pore pressure decreases in the case when the tunnel is placed in the trench.

Figure 9 shows the variation of the pore pressure at 1 m below the bottom two corners of the tunnel with time. The pore pressure amplitude on the lee side is shown to be larger than that on the weather side, which may be caused by the changes of water depth. While in Figure 9b, the amplitude of the pore pressure on the weather side is larger than that on the lee side. Compared with Figure 9a, the amplitude in Figure 9b on the lee side is reduced slightly. The violent change of pore pressure may lead to a more serious transient liquefaction.

3.4. Wave-Induced Liquefaction

The wave-induced liquefaction of soil is important for the design and construction of offshore engineering structures. Meanwhile, several wave-induced liquefaction criterions have been put forward [35,36,37]. Zen and Yamazaki [35] proposed a two-dimensional liquefaction criterion:

$${\sigma }_{v0}^{\prime }<p_s-p_b$$

(27)

in which $p_b$ is the wave pressure acting on the seabed surface and $p_s$ is the wave-induced oscillatory pore pressure in the seabed.

Figure 10 presents the liquefaction in the vicinity of the trench at the same time (t = 36 s) in two cases: (a) after the foundation groove is excavated and (b) after the tunnel is settled down. The maximum liquefied depth in case (b) is 0.95 m, smaller than that in case (a) with 1.75 m. In the wave field, the maximum wave pressure in case (b) is 1.90 × 10⁴ Pa, smaller than the 2.01 × 10⁴ Pa in case (a). This may be due to a certain resistance to the wave propagation from the tunnel.

3.5. Wave-Induced Shear Failure

Generally speaking, when the shear stress is greater than the shear strength of the soil, the shear failure occurs and the relative displacement of the soil particles leads to instability failure [38]. Although the Mohr–Coulomb criterion is based on the linear elastic theory, it still is a basic and convenient tool for engineering reference and is used here to determine the shear failure of the seabed. First, it is necessary to comprehend the effective stress path of the soil. The $p^{\prime }-q^{\prime }$ plane is adopted to analyze the effective stress path; $p^{\prime }$ and $q^{\prime }$ are respectively defined as

$$p^{\prime }=({\sigma }_1^{\prime }+{\sigma }_3^{\prime })/2$$

(28)

$$q^{\prime }=({\sigma }_1^{\prime }-{\sigma }_3^{\prime })/2$$

(29)

where ${\sigma }_1^{\prime }$ and ${\sigma }_3^{\prime }$ are the effective principal stresses in the seabed, respectively. The effect of intermediate principal stress is not considered.

Thus, considering the initial consolidation state of the seabed, the initial stress without a wave load can be expressed as

$${\sigma }_{z0}^{\prime }={\gamma }^{\prime }z$$

(30)

$${\sigma }_{x0}^{\prime }=K_0{\gamma }^{\prime }z$$

(31)

where $K_0$ is the coefficient of the lateral earth pressure, which is related to Poisson’s ratio; $K_0=\mu /(1-\mu )$; and ${\gamma }^{\prime }$ is the effective unit weight of soil.

Herein, after being subjected to the wave load, the total effective stress of soil can be expressed as

$${\overline{\sigma }}^{\prime }{}_z={\sigma }_{z0}^{\prime }+{\sigma }_{zz}^{\prime }$$

(32)

$${\overline{\sigma }}^{\prime }{}_x={\sigma }_{x0}^{\prime }+{\sigma }_{xz}^{\prime }$$

(33)

where ${\sigma }_{zz}^{\prime }$ and ${\sigma }_{xz}^{\prime }$ denote the vertical and horizontal effective stresses in soil under wave loads.

Since the shear stress under the initial dead weight is equal to 0, the final total shear stress ${\overline{\tau }}^{\prime }{}_{xz}$ is induced by the wave load, which is expressed as

$${\overline{\tau }}^{\prime }{}_{xz}={\tau }_{xz}^{\prime }$$

(34)

The effective principal stresses ${\sigma }_1^{\prime }$ and ${\sigma }_3^{\prime }$ are expressed as

$${\sigma }_1^{\prime }=\frac{{\overline{\sigma }}^{\prime }{}_x+{\overline{\sigma }}^{\prime }{}_z}{2}+\sqrt{{\left(\frac{{\overline{\sigma }}^{\prime }{}_x-{\overline{\sigma }}^{\prime }{}_z}{2}\right)}^2+{\left({\overline{\tau }}^{\prime }{}_{xz}\right)}^2}$$

(35)

$${\sigma }_3^{\prime }=\frac{{\overline{\sigma }}^{\prime }{}_x+{\overline{\sigma }}^{\prime }{}_z}{2}-\sqrt{{\left(\frac{{\overline{\sigma }}^{\prime }{}_x-{\overline{\sigma }}^{\prime }{}_z}{2}\right)}^2+{\left({\overline{\tau }}^{\prime }{}_{xz}\right)}^2}$$

(36)

Then, the stress statement at one point in the soil part can be expressed as

$$\sin \varphi =\frac{{\sigma }_1^{\prime }-{\sigma }_3^{\prime }}{{\sigma }_1^{\prime }+{\sigma }_3^{\prime }+2c/\tan {\varphi }_f}$$

(37)

where $\varphi$ is the stress angle, and $c$ and $\varphi$ are the cohesion and internal friction angles, respectively.

Based on the Mohr–Coulomb criterion, the shear strength of the soil can be expressed as

$${\tau }_f={\sigma }_f\tan {\varphi }_f+c$$

(38)

where ${\sigma }_f$ and ${\tau }_f$ are the normal stress and shear stress on the failure surface, separately.

When the stress in the soil reaches the strength envelope, the soil reaches the ultimate equilibrium state. Therefore, the discriminant formula of shear failure at one point in the soil can be written as

$$\varphi \ge {\varphi }_f$$

(39)

Substitute Equation (39) into Equation (37), and the discriminant formula of shear failure is obtained as

$$\mathrm{a}\varphi \ge \arcsin \left(\frac{\frac{{\sigma }_1^{\prime }-{\sigma }_3^{\prime }}{2}}{\frac{c}{\tan \varphi }+\frac{{\sigma }_1^{\prime }+{\sigma }_3^{\prime }}{2}}\right)\ge {\varphi }_f$$

(40)

For a more detailed introduction about shear failure, readers can refer to Zen et al. (1998) [39] and Jeng (2001) [38].

Figure 11 shows the distribution of the shear failure area in the trench in two cases: (a) after the foundation groove is excavated and (b) after the tunnel element is placed. The results at t = 36 s are presented as the depth of the shear failure area when it reaches the maximum at this moment. The maximum depths of shear failure at the left toe are almost the same (at z = −33.5 m) in these two cases. Due to the existence of the tunnel, the depth of the shear failure area below the tunnel reduces. Moreover, the shear failure is more likely to occur near two corners at the bottom of the tunnel because the stress is concentrated here.

3.6. Influence of Wave Characters on Liquefaction

In this section, based on the proposed model, the influences of the wave parameters on the dynamic response of soil around the foundation trench are studied. The standard wave conditions are $T=8 \mathrm{s}$, $H=4 \mathrm{m}$ and $d=16 \mathrm{m}$. In the parametric analysis, the wave periods are $T=6 \mathrm{s}$, $7 \mathrm{s}$ and $8 \mathrm{s}$; the wave heights are $H=3 \mathrm{m}$, $4 \mathrm{m}$ and $5 \mathrm{m}$; and the water depths are $d=12 \mathrm{m}$, $16 \mathrm{m}$ and $20 \mathrm{m}$. The calculation domain is 100 < x < 200 m and −46.5 < y < −16 m, and the maximum depth of liquefication, the area of liquefication zone and the width of liquefication at the seabed surface are calculated.

In general, the waves with a large wave length and a large wave height in shallow water tend to induce liquefaction in the seabed easily [6]. As shown in Figure 12, with the increases of wave height $H$, the increases of wave period $T$ or the decreases of water depth $d$, the liquefaction depth, the liquefaction area and the liquefaction width increase. However, the increase of the liquefaction area is mainly associated with the increase of the liquefied depth. Moreover, the increase of water depth can effectively restrain the dynamic pressure generated by waves, and thus, the dynamic response would be much smaller.

3.7. Influence of Slope Rate on Liquefaction

Then, the influence of the slope side angle to the stability of the foundation trench is further investigated. It is obvious that once the bottom width and the excavation depth of the trench are fixed, the factor affecting the amount of excavation and backfill is exactly the slope angle $\theta$. The suitable $\theta$ could ensure the stability of the excavated trench, while meeting the economic requirements. In this section, the wave adopts the standard condition from above and the slope ratios are taken to be 1:3, 1:2 and 1:1.5, respectively. As discussed previously, the maximum liquefaction depth, the liquefaction area and the liquefaction width are calculated, and the latter two are expressed in the forms of percentages (the liquefied area is nondimensionalized using the area of the trench and the liquefied width using the width of trench bottom). As shown in Figure 13, with the increase of the slope ratio, the maximum liquefaction depth, the liquefaction area and the liquefaction width decrease at first and then increase without exception. That is to say, under certain wave conditions, there is an optimal value for the slope ratio.

4. Conclusions

In this paper, a two-dimensional coupling model of a wave-seabed-immersed tunnel is proposed, which aims to study the dynamic response of the soil around the foundation trench and immersed tube tunnel under wave load. The model is verified using the analytical and experimental results. The influences of the wave characteristics, slope ratio and soil parameters on the soil responses are analyzed. The influence of soil failure on the stability of the tunnel is also discussed. Based on the numerical results, the main conclusions can be drawn as follows:

(1)

Due to the existence of the trench, the pore pressure amplitude on the weather side slope is significantly smaller than that on the lee side slope.

(2)

The maximum depth of liquefaction in the case after the tunnel element is placed is smaller than that after the foundation groove is excavated.

(3)

Due to the existence of the tunnel structure, the distribution of the flow field and pressure field change dramatically; thus, the dynamic responses and the failure area in the seabed change accordingly.

(4)

In the case of the specific wave and seabed parameters, the liquefaction characteristics in the trench have an obvious fold point with the change of slope rate. That means that there is an optimal slope rate to minimize the failure possibility of the slope. Moreover, the specific failure mode deserves further research.