Integrated Model for Wave-Induced Oscillatory and Residual Soil Response in a Poro-Elastic Seabed: Partially Dynamic Model

Abstract

The evaluation of wave-induced residual pore pressure in a porous seabed and associated seabed liquefaction is essential for designing marine infrastructure foundations. The strength and stiffness of the seabed could be weakened due to the build-up of pore pressures under cyclic wave action, further leading to residual liquefaction. Existing models for residual liquefaction are limited to the quasi-static uncoupled approaches, which do not account for the effect of oscillatory pore pressure on the accumulative pore pressure acceleration of solid particles, despite the mutual influence of these two mechanisms. To overcome these limitations, this paper proposes a new model for residual soil response with $u-p$ approximation (partial dynamic model) that couples oscillatory and residual mechanisms. The proposed model is validated through wave flume tests and centrifuge tests. Based on the coupling model, a new criterion of liquefaction integrating both oscillatory and residual mechanisms is also proposed. Numerical examples demonstrate that the coupling effect significantly affects the wave-induced seabed liquefaction potential. Furthermore, a new parameter ($\mathsf{\Omega }$) representing the ratio of oscillatory and residual pore pressure is introduced to clarify which mechanism dominates the pore pressure development.

1. Introduction

The pore pressure within a seabed is a crucial parameter that directly impacts the design of foundations of marine infrastructures. Evaluating soil characteristics and seabed stability by using pore pressure have practical engineering applications, especially in newly deposited seabeds with soft and in-compact soil, which are prone to instability and low bearing capacity [1,2]. Two wave-induced seabed liquefaction mechanisms, namely oscillatory and residual, have been reported in the literature [3,4,5]. The former is caused by oscillatory or transient excess pore pressures with phase lag and amplitude attenuation [6,7,8,9,10,11,12], while the latter is the result of the build-up of excess pore pressures due to soil contraction under cyclic loading [4,13]. This study considers both mechanisms.

Previous studies on wave-induced residual liquefaction were based on Biot’s consolidation model, a quasi-static approach that neglects accelerations caused by soil and fluid particles [13,14,15,16,17,18]. However, soil and fluid particle accelerations are critical under specific combinations of wave and seabed conditions [19,20], which are not considered in the decoupled approach used in previous studies. Based on Biot [21] and Zienkiewicz et al. [22], Jeng et al. [23] proposed the first $u-p$ approximation for the wave-induced oscillatory soil response. However, to date, there is no $u-p$ approximation for the wave-induced residual liquefaction available in the literature.

In addition to the acceleration of soil particles, the previous studies for the wave-induced residual pore pressure adopted decoupled approaches. In the approach, the oscillatory seabed response is calculated for the whole spatial and time domains, and then the shear stress of the oscillatory mechanism is employed to determine the source term for the generation of the residual pore pressures. This approach assumes that the oscillatory soil shear stress and source term do not change with the development of the residual pore pressure. The major limitation of the aforementioned models is that the oscillatory and residual components of the total pore pressure do not affect each other. That is, the coupling effect between residual and oscillatory mechanisms is ignored. In fact, it has been confirmed that the time series of the pore water pressure includes both oscillatory and residual components [4], and soil shear stress and pore pressure are related under wave loading [24,25]. Recently, Liu et al. [26] explored the difference between coupled and uncoupled approaches regarding the pore pressure based on a quasi-static model.

This paper will make two new contributions to the field of the wave-seabed interactions. They are

• First, a new model for wave-induced residual pore pressure based on $u-p$ approximation will be proposed, in which the acceleration of soil particles will be included in the model for residual mechanism.
• Second, the oscillatory and residual mechanisms will be coupled, in which both mechanisms will affect each other. Furthermore, this study will clarify when individual mechanisms will dominate the development of pore pressure.

The paper is organised as follows: Section 2 outlines the theoretical model, including the oscillatory and residual seabed models and the coupling process, which is established within OpenFOAM. Section 3 validates the present model by comparing it with previous experimental data and the decoupled model [18]. Section 4 presents a detailed discussion of both mechanisms based on the present model. Finally, Section 5 summarizes the key findings.

2. Theoretical Model

Figure 1 illustrates the interactions between propagating ocean waves and a porous seabed. As shown in the figure, the waves travel in the positive x-axis, while the z-axis is upward from the seabed surface. In this study, the linear wave theory is adopted as the first approximation [27], in which the wave profile ($\eta$) and dynamic wave pressure ($p_b$) can be expressed as:

$$\eta =\frac{H}{2}cos(kx-\omega t),$$

(1)

$$p_b=\frac{{\gamma }_{w}H}{2}\frac{coshkz}{coshkd}cos(kx-\omega t),$$

(2)

where $\eta$ is the free water surface, $p_b$ is the dynamic wave pressure, ${\gamma }_w$ is the unit weight of water, $k=2\pi /L$ is the wave number, $\omega =2\pi /T$ is the wave frequency (T represents the wave period), H is wave height, and d is wave depth.

Based on the wave dispersion relation, the wavelength (L) can be obtained using

$$L=\frac{g{T}^2}{2\pi }tanh\left(\frac{2\pi d}{L}\right).$$

(3)

The present seabed model is established under the OpenFOAM platform. The partially dynamic model ($u-p$ approximation) is adopted for both oscillatory and residual soil response for the development of the seabed model. The seabed is assumed to be hydraulically isotropic with the same permeability ($k_s$) in all directions, in which the pore fluid is considered compressible and obeys Darcy’s law.

Figure 2 illustrates the development of wave-induced pore pressures and two mechanisms, namely oscillatory (${\tilde{p}}_s$) and residual (${\overline{p}}_s$). The pore pressure ($p_s$) can be expressed as [4,28]

$$p_s={\tilde{p}}_s+{\overline{p}}_s,{\overline{p}}_s=\frac{1}{T}{\int }_{t_0}^{t_0+T}{p}_{s}dt.$$

(4)

2.1. Partially Dynamic ($u-p$) Oscillatory Seabed Model

Considering the oscillatory mechanism in two-dimension (2D), the equation that governs the force equilibrium in soil is shown as [23]:

$$G_s{\nabla }^2{\mathbf{u}}_s+\frac{G_s}{1-2{\mu }_s}\nabla (\nabla ·{\mathbf{u}}_s)=\nabla p_s+\rho \frac{{\partial }^2{\mathbf{u}}_s}{\partial t^2},$$

(5)

where ${\nabla }^2$ represents the Laplace operator; ${\mathbf{u}}_s\equiv (u_s,w_s)$ is the soil displacement vector; $G_s$ is shear modulus, ${\mu }_s$ is the Poisson’s ratio, $p_s$ is pore pressure, and $\rho$ is the mixture density of soil and water, expressed as $\rho =n_s{\rho }_w+(1-n_s){\rho }_s$, where ${\rho }_w$ is water density and ${\rho }_s$ is soil density. Note that the last term on the right-hand side represents the acceleration of soil particles, which only appears in the partial dynamic model.

The effective stress-strain relation is based on the generalised Hooke’s law:

$${\sigma }_x^{\prime }=2G_s\left[\frac{\partial u_s}{\partial x}+\frac{{\mu }_s(\nabla ·{\mathbf{u}}_s)}{1-2{\mu }_s}\right],$$

(6)

$${\sigma }_z^{\prime }=2G_s\left[\frac{\partial w_s}{\partial z}+\frac{{\mu }_s(\nabla ·{\mathbf{u}}_s)}{1-2{\mu }_s}\right],$$

(7)

$${\tau }_{xz}=G_s\left[\frac{\partial u_s}{\partial z}+\frac{\partial w_s}{\partial x}\right]={\tau }_{zx},$$

(8)

in which ${\sigma }_x^{\prime }$ and ${\sigma }_z^{\prime }$ are effective normal stresses and ${\tau }_{xz}$ is the shear stress. $G_s$ is the shear modulus of soil: $G_s=E_s/2(1+{\mu }_s)$, and $E_s$ is the Young’s modulus.

Based on the conservation of mass, the governing equation for the pore fluid,i.e., storage equation, can be expressed as [23]:

$${\nabla }^2{p}_s-\frac{{\gamma }_w{n}_s{\beta }_s}{k_s}\frac{\partial p_s}{\partial t}+\rho \frac{{\partial }^2(\nabla ·{\mathbf{u}}_s)}{\partial t^2}=\frac{{\gamma }_w}{k_s}\frac{\partial (\nabla ·{\mathbf{u}}_s)}{\partial t},$$

(9)

where $n_s$ and $k_s$ are soil porosity and permeability, respectively. Note that the last term on the left-hand-side of (9) only appears in the $u-p$ approximation.

In (9), ${\beta }_s$ is the pore fluid compressibility, which is defined as:

$${\beta }_s=\frac{1}{K_w}+\frac{1-S_r}{{\gamma }_{w}d}$$

(10)

where $K_w$ is the apparent bulk modulus of pore water, usually taking the value of 1.95 × 10${}^9$ N/m${}^2$ [6].

Equations (5) and (9) consist of the governing equations of the 2D partially dynamic ($u-p$) oscillatory seabed model. Based on the oscillatory model, the wave-induced oscillatory pore water pressure (${\tilde{p}}_s$), the oscillatory volume strain ($\widetilde{ϵ}=(\nabla ·{\mathbf{u}}_s)$) can be obtained.

In the following section, to distinguish the difference between oscillatory and residual soil responses, variables regarding the oscillatory part of pore pressure will be denoted with ’∼’.

2.2. Partially Dynamic ($u-p$) Residual Seabed Model

Considering the tension as positive and compression as negative for normal stresses (i.e., $\sigma ={\sigma }^{\prime }-p_s$), total normal stresses in horizontal and vertical directions, respectively, are shown as follows:

$${\sigma }_x={\sigma }_x^{\prime }-p_s=2G_s\left[\frac{\partial {\tilde{u}}_s}{\partial x}+\frac{{\mu }_s}{1-2{\mu }_s}(\nabla ·{\tilde{\mathbf{u}}}_s)\right]-p_s,$$

(11)

$${\sigma }_z={\sigma }_z^{\prime }-p_s=2G_s\left[\frac{\partial {\tilde{w}}_s}{\partial z}+\frac{{\mu }_s}{1-2{\mu }_s}(\nabla ·{\tilde{\mathbf{u}}}_s)\right]-p_s,$$

(12)

where $p_s$ is the pore pressure that consists of oscillatory component (${\tilde{p}}_s$) and residual component (${\overline{p}}_s$). In most existing models, these two mechanisms were considered individually. In this study, these two mechanisms will be integrated.

The average of the total normal stresses in the $x-$ and $z-$ direction can be expressed as

$${\sigma }_m=\frac{{\sigma }_x+{\sigma }_z}{2}=\frac{G_s}{1-2{\mu }_s}(\nabla ·{\tilde{\mathbf{u}}}_s)-p_s=\frac{1}{m_v}(\nabla ·{\tilde{\mathbf{u}}}_s)-p_s,$$

(13)

where $m_v$ is the volume compressibility coefficient:

$$m_v=\frac{1-2{\mu }_s}{G_s}$$

(14)

Based on Seed and Rahman [13] and Tsotsos et al. [29], change of the volumetric strains ($\mathsf{\Delta }ϵ=\mathsf{\Delta }(\nabla ·{\mathbf{u}}_s)$) in two dimensions can be expressed as

$$\mathsf{\Delta }ϵ=\mathsf{\Delta }({ϵ}_x+{ϵ}_z)=m_v(\mathsf{\Delta }{p}_s-\psi \mathsf{\Delta }t+\mathsf{\Delta }{\sigma }_m),$$

(15)

or

$$\frac{\partial ϵ}{\partial t}=m_v\left(\frac{\partial p_s}{\partial t}-\psi +\frac{\partial {\sigma }_m}{\partial t}\right)=m_v\left(\frac{\partial (p_s+{\sigma }_m)}{\partial t}\right)-m_v\psi ,$$

(16)

where $\mathsf{\Delta }{p}_s$ is the net pore pressure change during $\mathsf{\Delta }t$ calculated as the pore pressure generation minus pore pressure dissipation, and $\psi$ is a source term that represents the pore pressure generation rate under cyclic wave loading for undrained conditions. Here, we also consider the effect of mean total normal stress on pore pressure change, which is reflected in term $\mathsf{\Delta }{\sigma }_m$.

Combining (13) and (16), we have

$$\frac{\partial ϵ}{\partial t}=m_v\left[\frac{\partial p_s}{\partial t}-\psi +\frac{\partial }{\partial t}\left(\frac{1}{m_v}\widetilde{ϵ}-p_s\right)\right]=\frac{\partial \widetilde{ϵ}}{\partial t}-m_v\psi .$$

(17)

Substituting (17) into (9) yields the following:

$${\nabla }^2{p}_s-\frac{{\gamma }_w{n}_s{\beta }_s}{k_s}\frac{\partial p_s}{\partial t}+\rho \frac{{\partial }^2(\nabla ·{\tilde{\mathbf{u}}}_s)}{\partial t^2}=\frac{{\gamma }_w}{k_s}\frac{\partial (\nabla ·{\tilde{\mathbf{u}}}_s)}{\partial t}-\frac{{\gamma }_w}{k_s}{m}_v\psi ,$$

(18)

in which $\psi$ can be determined by [13];

$$\psi =\frac{\partial u_g}{\partial t}=\frac{{\sigma }_0^{\prime }}{T}{\left(\frac{\left|{\tau }_{max}\right|}{{\alpha }_r{\sigma }_0^{\prime }}\right)}^{-1/{\beta }_r},$$

(19)

where $u_g$ is the generation of pore pressure, ${\tau }_{max}$ is the maximum shear stress amplitude, and ${\sigma }_0^{\prime }$ is the mean value of initial effective stress, which is defined as

$${\sigma }_0^{\prime }=\frac{1+2k_0}{3}{\gamma }^{\prime }\left|z\right|,$$

(20)

where $k_0$ is lateral earth pressure coefficient; and ${\gamma }^{\prime }$ is the submerged unit weight of the soil.

In (19), ${\alpha }_r$ and ${\beta }_r$ are two residual parameters that are functions of the relative density of soil, which can be determined from the following empirical expressions [15,30]:

$${\alpha }_r=0.34D_r+0.084,{\beta }_r=0.37D_r-0.46,$$

(21)

where $D_r$ is the soil relative density, which can be determined by:

$$D_r=\frac{e_{max}-e}{e_{max}-e_{min}},$$

(22)

where e is the void ratio, and $e_{max}$ and $e_{min}$ are the maximum and minimum void ratios, respectively.

Note that the third term on the left-hand-side of (18) is an additional term comes from the $u-p$ approximation, which was not included in the previous residual model [4,18,31].

2.3. Boundary Conditions

Equations (5) and (9) together with (18) are the governing equations of the partially dynamic ($u-p$) residual mechanism. With appropriate boundary conditions, the wave-induced soil response can be obtained by solving the above governing equations with appropriate boundary conditions.

At the surface of the seabed ($z=0$), the pore pressure ($p_s$) is equal to the value of hydrodynamic pressure ($p_b$) obtained from the wave model, (i.e., (2)), and the effective normal stress and shear stresses are considered zero:

$$p_s=p_b(x,0,t)=\frac{{\gamma }_{w}H}{2coshkd}cos(kx-\omega t)=P_{0}cos(kx-\omega t),{\sigma }_z^{\prime }={\tau }_{xz}=0,atz=0,$$

(23)

where $P_0$ represents the amplitude of dynamic wave pressure at the seabed surface.

The seabed bottom ($z=-h$) is considered as impermeable rigid boundary, in other words, the vertical displacement and vertical flow gradient are zero:

$${\mathbf{u}}_s=\frac{\partial p_s}{\partial z}=0,\mathrm{at}z=-h$$

(24)

For the lateral boundaries, zero flow in the horizontal direction is applied. A large computational domain with a length that is over three times the wavelength is adopted in order to eliminate the impact from wave reflection. Such size of the computational domain has been proved to be sufficient for the seabed simulation [32].

2.4. Integrated Process of Two Mechanisms

To integrate both mechanisms (i.e., oscillatory and residual), the following steps will be taken.

(1)

In the first wave cycle, the oscillatory soil response (${\tilde{p}}_s^{\left(1\right)}$, etc) will be determined from (5)–(9) with appropriate boundary conditions at the first wave cycle. The maximum shear stress in the first wave cycle can be determined, which will be used to determine the source term using (19).

(2)

The residual pore pressure for the first wave cycle (${\overline{p}}_s^{\left(1\right)}$) is determined by solving (18).

(3)

In the next wave cycle, substituting $p_s^{\left(2\right)}={\tilde{p}}_s^{\left(2\right)}+{\overline{p}}_s^{\left(1\right)}$, into (5)–(9), we can determine the oscillatory soil response in the second wave cycle (${\tilde{p}}_s^{\left(2\right)}$), and obtain the maximum shear stress and source term for the second wave cycle. Then, the residual pore pressure for the second wave cycle (${\overline{p}}_s^{\left(2\right)}$) can be determined using (18). This will integrate both oscillatory and residual mechanisms.

(4)

Repeating the above procedure for the following wave cycles by replacing $p_s^{\left(1\right)}$ and $p_s^{\left(2\right)}$ with $p_s^{\left(n\right)}$ and $p_s^{(n+1)}$.

3. Mesh Convergence

Figure 3 and Figure 4 show the mesh convergence analysis for mesh size in the horizontal and vertical directions, respectively. In the figures, a parameter indicating the computational efficiency ($CE$) is also included, which is defined as the actual execution time divided by 100 s. That is, $CE=$ 400 means the execution time for the mesh is 40,000 s (CPU time). The pore pressure calculated using the finest mesh ($\mathsf{\Delta }x$ = 0.354, $\mathsf{\Delta }z$ = 0.1) is adopted as the reference value to test the convergence of other meshes, which are named $p_{ref}$ in figures. As it can be seen from the figures, when $\mathsf{\Delta }x$ is 0.354 and $\mathsf{\Delta }z$ is 0.1, the simulation can be considered convergent with the best computational efficiency, as marked in red in the figures. Therefore, the grid of the numerical model in the following numerical examples will adopt the above mesh scheme.

4. Model Validation

In order to validate the present model, two sets of experimental data are adopted. The first experimental data is the wave flume tests conducted by [15], while the second experimental data is the centrifuge test using [33]. In addition to the comparisons with experimental data, the quasi-static (QS) de-coupled model of residual soil response that developed within OpenFOAM [18] is also included in the comparisons.

4.1. Validation #1: Comparison with Wave Flume Tests [15]

Sumer et al. [15] carried out wave flume tests to investigate the development of pore pressure in a silty seabed and validate the accuracy of the quasi-static decoupled model. The test settings included a piston-type wave generator, the wave absorber, and sedimentary soil pit. The water depth was always maintained at 55 cm, while the pit was 14 m from the wave generator. Pore pressure transducers and wave gauges were set up to monitor the pore pressure and wave surface elevations, respectively. The parameters are identical to the wave flume tests, which are listed in

Table 1. Parameters for the first validation [15].

Characteristics	Value	Unit
Wave characteristics
Wave height (H)	0.18	[m]
Water depth (d)	0.55	[m]
Wave period (T)	1.6	[s]
Wave length (L)	3.18	[m]
Bulk modules (K)	1.9 × 10${}^9$	[N/m${}^2$]
Amplitude of pore pressure ($P_0$)	534.63	[pa]
Seabed characteristics
Porosity ($n_s$)	0.51	–
Permeability ($k_s$)	1.5 × 10${}^{-5}$	[m/s]
Seabed length ($L_s$)	3L	[m]
Poisson’s ratio (${\mu }_s$)	0.29	–
Shear modules ($G_s$)	1.92 × 10${}^6$	[N/m${}^2$]
Relative density ($D_r$)	0.28	–
Seabed thickness (h)	0.4	[m]
Degree of saturation ($S_r$)	100	%
Submerged weight of soil (${\gamma }^{\prime }$)	8140	[N/m${}^3$]
Pressure buildup characteristics
Residual parameter (${\alpha }_r$)	0.1792	–
Residual parameter (${\beta }_r$)	−0.3564	–
Lateral earth pressure coefficient ($k_0$)	0.42	–
Coefficient of consolidation ($c_v$)	9.71 × 10${}^{-3}$	–

.

Figure 5 shows that the results from the present integrated model are in good agreement with wave flume test data. However, the pore pressure simulated from the decoupled model [18] is not fully developed at the initial stage of loading, especially in Figure 5b. The present model has the ability to reproduce both oscillatory and residual components of the pore pressure.

The main differences between the test and the present model are the phase lag and the increasing speed of pore pressure. There are two reasons for this phenomenon: one is the error of the experimental measurements, and the other is the selection of the source term, which account for the difference between the experimental data and the models in Figure 5b more than in Figure 5a. In the process of pore pressure accumulation, especially at the initial stage of loading, i.e., 0∼5 s, it is possible to cause the pore pressure accumulation to be faster due to the determination of the source term corresponding to the maximum shear stress in the oscillatory model. In order to further illustrate the fitting accuracy, RMSE (root mean square error) is used to analyse differences between simulations and tests, which is defined as

$$RMSE=\sqrt{\frac{1}{N}\sum _{i=1}^N{(p_{num}-p_{exp})}^2},$$

(25)

where $p_{num}$ is the numerical results; $p_{exp}$ is the experimental data and N is the number of data points.

Based on Figure 5, 50 experimental data points are selected in wave loading periods, and the simulated pore pressure values corresponding to the same time in the tests are taken in. Finally, the RMSE values are listed in

Table 2. Comparison of the RMSE of the present coupling model and decoupled model [18] with experimental data [15].

Location	Coupling Model	Decoupled Model
$(x,z)$ = (5.0, −0.085) m [Figure 5a]	1.32	4.14
$(x,z)$ = (5.0, −0.24) m [Figure 5b]	5.38	7.56

. As shown in the table, the u-p coupling model can provide better accuracy for simulating real experiments, compared with decoupled model [18].

4.2. Validation #2: Comparison with Centrifugal Test [33]

The second validation case is to compare with data from Sassa and Sekiguchi [33]’s centrifuge wave tests, which was to investigate the soil liquefaction under standing and progressive waves. The centrifuge tests were carried out under 50× g acceleration, in the numerical set-up, the parameters with a scaling factor of 50 are used as suggested by Sassa and Sekiguchi [33]. The parameters are listed in

Table 3. Parameters for the second validation [18].

Characteristics	Value	Unit
Wave characteristics
Wave height (H)	2.06	[m]
Water depth (d)	4.5	[m]
Wave period (T)	4.5	[s]
Wave length (L)	25.43	[m]
Bulk modules (K)	1.9 × 10${}^9$	[N/m${}^2$]
Amplitude of pore pressure ($P_0$)	6000	[pa]
Seabed characteristics
Porosity ($n_s$)	0.45	–
Permeability ($k_s$)	1.5 × 10${}^{-4}$	[m/s]
Seabed length ($L_s$)	3L	[m]
Poisson’s ratio (${\mu }_s$)	0.3	–
Shear modules ($G_s$)	1.0 × 10${}^6$	[N/m${}^2$]
Relative density ($D_r$)	0.437	–
Seabed thickness (h)	5.0	[m]
Degree of saturation ($S_r$)	100	%
Submerged weight of soil (${\gamma }^{\prime }$)	8600	[N/m${}^3$]
Pressure buildup characteristics
Residual parameter (${\alpha }_r$)	0.2326	–
Residual parameter (${\beta }_r$)	−0.2983	–
Lateral earth pressure coefficient ($k_0$)	0.52	–
Coefficient of consolidation ($c_v$)	5.25 × 10${}^{-2}$	–

.

Figure 6 presents the pore pressure build-up from the present model, centrifugal test, and decoupled model [18]. As can be seen in the figure, at the depth of 2.5 m, the result from the present model is in good agreement with experimental data. Similar to the last case, decoupled model underestimates the build-up of pore pressure. Compared with the decoupled model, the present model can predict more accurately the pore pressure build-up in the centrifugal tests.

5. Results and Discussions

5.1. Cases Considered in Numerical Examples

In this section, the development of the wave-induced pore pressure and the associated liquefaction will be discussed. The parameters selected in the following numerical examples can be found in

Table 4. Parameters used in numerical examples.

Case Number	H (m)	d (m)	T (s)	$\mathit{H}/\mathit{L}$	$\mathit{d}/\mathit{L}$	${\mathit{k}}_{\mathit{s}}$ (m/s)
1	2.8	4.0	5.0	0.100	0.143	1.0 × 10${}^{-4}$
2	2.8	4.0	4.0	0.134	0.191	1.0 × 10${}^{-4}$
3	2.8	4.0	4.0	0.134	0.191	7.0 × 10${}^{-4}$

Note: other parameters, $n_s$ = 0.45; ${\mu }_s$ = 0.35; $G_s$ = 3.0 × 10${}^6$; $D_r$ = 0.54; h = 20 m; $k_0$ = 0.42.

. As reported in the literature, an oscillatory mechanism is likely to dominate the development of pore pressure in Case 1, while Case 2 will be dominated by a residual mechanism based on the de-coupled model [14]. When the pore pressure coupling model is dominated by an oscillatory mechanism in the first case, the progressive wave loading period is 5.0 s, the wavelength is 27.9 m, and the amplitude of pore pressure is 9.59 kPa for Case 1, while the wave period is 4.0 s, the wavelength is 20.9 m and the amplitude of pore pressure is 7.55 kPa for Case 2. The parameters used in the third case (Case 3) are consistent with Case 2, except for the different permeability. The middle section of the seabed (at $x=50$ m) is used as a monitoring position for analysing the pore pressure development.

5.2. Liquefaction Criterion

The wave-induced pore pressure has been used in the evaluation of the liquefaction in a porous seabed [4]. Herein, we outline the criteria for momentary and residual liquefaction first and then propose a new criterion to include both mechanisms.

The momentary liquefaction is caused by oscillatory soil response, which can be determined using the following criterion [3]:

$${\tilde{p}}_s-p_b\ge {\sigma }_0^{\prime }.$$

(26)

One the other hand, the residual liquefaction can be determined using the following criterion [4,16]:

$${\overline{p}}_s\ge {\sigma }_0^{\prime }.$$

(27)

In the existing studies for the wave-induced liquefaction in a porous seabed, either (26) or (27) is used in the liquefaction assessment. That is, momentary liquefaction and residual liquefaction are considered separately. In this study, we integrate both oscillatory and residual mechanisms by considering that the pore pressure fluctuates while it is building-up. Therefore, the above criteria of liquefaction need to merge and form the new criterion of liquefaction to consist of both mechanisms, that is,

$$({\tilde{p}}_s-p_b)+{\overline{p}}_s=p_s-p_b\ge {\sigma }_0^{\prime }.$$

(28)

Based on the new criterion of liquefaction, together with the input data in Table 4, the following numerical examples will be presented in the discussions. It is observed that the nominal pore pressure ($(p_s-p_b)/P_0$, where $P_0$ is defined in (23)) is developing gradually with the wave cycles ($t/T$) both in Figure 7 and Figure 8. As it can be easily observed in Figure 7, in Case 1, compared with the slight accumulation of pore pressure over 30 wave cycles (i.e., the pore pressure accumulation is around 0.4 $P_0$ for location at a 2-m seabed depth), the amplitude of oscillating pore pressure accounts for a much larger proportion (i.e., the oscillating pore pressure amplitude is around 1.2 $P_0$ for location at a 2-m seabed depth), in other words, Case 1 is dominated by oscillatory mechanism. On the other hand, as shown in Figure 8 for Case 2, the amount of cumulative pore pressure is around 12 $P_0$ at a 2-m seabed depth over 70 wave cycles, while the amplitude of oscillating pore pressure is less than 2 $P_0$ at the same location. Therefore, under the condition of wave loading in Case 2, the same seabed will be dominated by residual mechanism rather than oscillatory mechanism. It is interesting to notice that the change in wave period could have a significant impact on the pore pressure response in the seabed with the same soil properties.

In both figures, the values of ${\sigma }_0^{\prime }/P_0$ at different locations are also presented for the purpose of indicating the liquefaction state. When ($(p_s-p_b)/P_0$) is greater than ${\sigma }_0^{\prime }/P_0$, liquefaction occurs. In Case 1 (Figure 7), after about 30 wave cycles, the cumulative pore pressure will reach the maximum where the liquefaction range no longer increases. However, in Case 2 (Figure 8), the cumulative pore pressure develops rapidly with time, which far exceeds the values of ${\sigma }_0^{\prime }/P_0$. It can be further seen that the seabed liquefies at z = −0.5 and −1.0 m (i.e., Figure 7a,b) but remains un-liquefied at z = −2.0 and −3.0 m (i.e., Figure 7c,d) in Case 1, which oscillatory liquefaction that appears under wave trough and disappears under wave crest is the dominance. While the seabed liquefies at all four locations in Case 2, the residual liquefaction is the dominance.

The pore pressure development for Case 3 is shown in Figure 9, which has the same wave period as Case 2 but a different soil permeability of 7.0 × 10${}^{-4}$ m/s. It is a case somewhere between Case 1 and Case 2, in which both the oscillatory mechanism and residual mechanism play an important role. For example, for location near the seabed surface (i.e., z = −0.5 m and −1.0 m), the residual pore pressure plays a leading role for seabed liquefaction, in other words, the seabed liquefies regardless of oscillatory pore pressure. While for the deeper locations, especially at z = −3.0 m, two mechanisms are equally important. Considering only either one of them solely may cause a miscalculation of the liquefaction.

Figure 10 and Figure 11 further illustrate the development of the liquefaction zone of the oscillatory-dominant and cumulative-dominant proposed model along the depth direction, respectively. As shown in Figure 10a–c, liquefaction mainly occurs from the equilibrium position to the trough position of the progressive wave. It can be observed in Figure 11 that liquefaction occurs mainly at position $z/h\ge -0.12$ in the second wave cycle. After about 10–30 cycles, the liquefied area develops to a certain degree at position $-0.35\le z/h\le -0.25$ in Figure 11b,c. Affected by the development of cumulative pore pressure in Figure 12a–c, the trend of the coupling liquefaction is similar to the cumulative-dominant cases with the relatively shallow liquefaction area.

To further compare the cases, Figure 13 illustrates the time series of the maximum liquefaction depth for Case 1, Case 2 and Case 3, respectively. As shown in the figure, the residual-dominant case (Case 2) has much larger maximum liquefaction depths and a faster accumulation rate compared to the oscillatory-dominant case, for example, at the 10th wave cycle, the maximum liquefaction depth reaches almost 5 m for Case 2 (residual-dominant) while only around 0.5 m for Case 1 (oscillatory-dominant). It indicates that a shorter wave period can achieve a deeper liquefaction area, resulting in residual-dominant liquefaction failure. In Case 3, the maximum liquefaction depth is always between that of Case 1 and Case 2, which is affected by both cumulative mechanism and oscillatory mechanism.

In Figure 13, the liquefaction depth of Case 2 is larger than Case 3. This is because Case 2 is dominated by the residual mechanism, in which the pore pressure reaches the initial effective stress very quickly, and the liquefaction depth grows along with time, whereas the residual mechanism for Case 3 has a lesser proportion compared with Case 2. The liquefaction criterion is reached based on the combination of the oscillatory and residual pore pressures, that is, the crest value of oscillatory on top of the residual pore pressures.

5.3. Momentary Liquefaction and Residual Liquefaction

To distinguish the influence of the development of cumulative components on the total pore pressure, the concept of cumulative share ratio is proposed to characterise the ratio of cumulative components to oscillatory components when the total pore pressure reaches a certain stable value.

The accumulative ratio that represents the ratio of residual and oscillatory pore pressure can be expressed as:

$$\mathsf{\Omega }=\frac{\left|{\overline{p}}_s\right|}{\left|{\tilde{p}}_s\right|}.$$

(29)

Note that the oscillatory mechanism will dominate the development of pore pressure when $\mathsf{\Omega }<<1$, while the residual mechanism will be more important when $\mathsf{\Omega }>>1$. Both mechanisms are equally important when $\mathsf{\Omega }=1$.

When $t/T=50$ is considered in three cases, the $\mathsf{\Omega }$ values are listed here.

$$\begin{array}{cc}\hfill Case1:& \mathsf{\Omega }=0.71,\mathrm{at}(x,z)=(50.0,-1.0)\mathrm{m}\hfill \\ \hfill Case2:& \mathsf{\Omega }=4.55,\mathrm{at}(x,z)=(31.0,-1.0)\mathrm{m}\hfill \\ \hfill Case3:& \mathrm{\Omega }=1.02,\mathrm{at}(x,z)=(31.0,-1.0)\mathrm{m}\hfill \end{array}$$

(30)

6. Conclusions

In this study, a $u-p$ integrated model is utilised to investigate the pore pressure accumulation in a porous seabed. Based on the simulation results, the following conclusions can be drawn:

(1)

In the case of oscillatory-dominant soil response (Case 1), the development of pore pressure is greatly affected by the oscillatory part of pore pressure. The maximum value of total pore pressure appears after 30 cycles of wave loading. Compared with the residual-dominant cases, the pore pressure growth rate is relatively slow. Resulting from the residual-dominant cases (Case 2), the development of pore pressure is greatly affected by residual pore pressure but less by oscillatory pore pressure. The pore pressure reaches the peak after 70 cyclic loading cycles, and the pore pressure increases rapidly.

(2)

The maximum liquefaction depth of the oscillatory-dominant liquefaction zone (Case 1) is relatively stable and will not change significantly with the wave cycle, while the development of the residual-dominant liquefaction depth (Case 2) is rapid, resulting in a large liquefaction range. The initial loading process of pore pressure cumulative liquefaction under the coupling mechanism (Case 3) is similar to that under the oscillatory mechanism. The maximum liquefaction trend of the coupling mechanism is similar to that under the residual mechanism, and the maximum liquefaction depth is between the two mechanisms.

(3)

Using the concept of accumulative ratio ($\mathsf{\Omega }$), the development trend of accumulative pore pressure is clarified. When $\mathsf{\Omega }$ approaches 1.0, residual pore pressure and oscillatory pore pressure play an equally important role. It may be an important feature for analysing the applications of coupling or decoupled models.

In this paper, only some preliminary results for the present coupling model are presented, and no structure is considered. In the future, the present model can be further applied to the evaluation of the wave-induced seabed liquefaction around marine infrastructure.