**Coastal Geohazard and Offshore Geotechnics**

Editors

**Dong-Sheng Jeng Jisheng Zhang V.S. Ozgur Kirca**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Dong-Sheng Jeng Griffith University Gold Coast Campus Australia

Jisheng Zhang Hohai University China

V.S. Ozgur Kirca Istanbul Technical University Turkey

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Journal of Marine Science and Engineering* (ISSN 2077-1312) (available at: https://www.mdpi.com/ journal/jmse/special issues/coastal geohazard).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-0274-8 (Hbk) ISBN 978-3-0365-0275-5 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**



## **About the Editors**

**Dong-Sheng Jeng** PhD, is currently a Professor at Griffith University, Gold Coast Campus, Australia. His research interest and expertise include theoretical modelling of fluid–seabed–structure interactions, porous flow, groundwater hydraulics, ocean/coastal engineering, offshore wind energy and artificial neural network modelling.

**Jisheng Zhang** PhD, received his PhD degree from the University of Aberdeen, Scotland, and is currently a Professor in Hohai University, College of Harbor, Coastal and Offshore Engineering, China. He has devoted himself to the research and application on seabed–flow–structure interactions, coastal dynamics and sediment, marine renewable energy engineering and environmental impacts.

**V.S. Ozgur Kirca** PhD, is currently an Associate Professor at Istanbul Technical University, Department of Civil Engineering, Turkey. He is also the co-founder of BM SUMER Consultancy and Research. His research interests and expertise include seabed–flow–structure interactions, turbulent processes, sediment transport and morphodynamics, physical and numerical modelling of coastal/hydraulic processes, and the design and assessment of coastal protection structures.

### *Editorial* **Coastal Geohazard and Offshore Geotechnics**

**Dong-Sheng Jeng 1,\*,†, Jisheng Zhang 2,† and Özgür Kirca 3,4,†**


Received: 5 December 2020; Accepted: 7 December 2020; Published: 10 December 2020

#### **1. Introduction**

With the rapid development in the exploration of marine resources, coastal geohazard and offshore geotechnics have attracted a great deal of attention from coastal geotechnical engineers and has achieved significant progress in recent years. With the complicated marine environment, numerous natural marine geohazard have been reported in the world, e.g., South China Sea. In addition, damage of offshore infrastructures (monopile, bridge piers, etc.) and supporting installations (pipelines, power transmission cables, etc.) have occurred in the last decades. A better understanding of the fundamental mechanism and soil behavior of the seabed in the marine environments will help engineers in the design or planning of the coastal geotechnical engineering projects. The purpose of this Special Issue is to present with the recent advances in the field of coastal geohazard and offshore geotechnics. This Special Issue will provide researchers updated development in the field and possible further developments.

In this Special Issue, eighteen papers were published, covering three main themes: (1) mechanism of fluid–seabed interactions and its associate seabed instability under dynamic loading [1–5]; (2) evaluation of stability of marine infrastructures, including pipelines [6–8], piled foundation and bridge piers [9–12], submarine tunnel [13], and other supported foundations [14]; and (3) coastal geohazard, including submarine landslide and slope stability [15,16] and other geohazard issue [17,18]. More details of each contribution are summarized in the following subsections.

#### **2. Mechanism and Processes of Seabed Response under Dynamic Loading**

The phenomenon of fluid–seabed interactions has attracted attentions among coastal and geotechnical engineers involved in the offshore geotechnical projects. A better understanding of the phenomenon and its associate processes will help practitioners and engineers in the design stage. The pore-water pressures and associated seabed liquefaction are key factors for the design of the foundation of offshore structures. The first theme of this Special Issue consists of five papers for the mechanism and processes of fluid–seabed interactions.

Liao et al. [1] proposed a coupling model for wave (current)-induced pore pressures and soil liquefaction in offshore deposits, based on the COMSOL Multiphysics. Unlike previous studies, both wave model and elastoplastic seabed model were established within COMSOL and coupled together, rather than through the data transformation at the fluid–seabed interface as the previous models. The numerical examples demonstrated the difference of the liquefaction depth between decoupling and coupling models. An alternative approach was proposed by Tong et al. [4] who integrated the

commercial software FLOW-3D and COMSOL Multiphysics for a similar problem, but with strong non-linear wave impact and uniform currents. More detailed discussions about the impact of current on the seabed response were provided.

Silty sand is a kind of typical marine sediment widely distributed in the offshore areas of East China. Guo et al. [3] investigated the wave-induced soil erosion in a silty sand seabed through a three-phase soil model (soil skeleton, pore fluid, and fluidized soil particles) within COMSOL Multiphysics. Based on their parametric study, it was found that the wave-induced erosion mainly occurred at the shallow depth of the seabed. Their study also found that the critical concentration of the fluidized soil particles has an obvious effect on the evolution of wave-induced erosion, including erosion rate and erosion degree. However, the erosion depth of seabed is not affected by the critical concentration of the fluidized soil particles.

Li et al. [5] integrated the hydrodynamic model (developed by OpenFOAM) and seabed model (developed by FEM) to investigate the effects of principal stress rotation (PSR) on the wave(current)-induced seabed liquefaction. The hydrodynamic model describes the process of the wave–current interactions. Meanwhile, the seabed model was based on the modified elastoplastic model with principal stresses. Based on their parametric study, it was found that PSR has significant impact on the development of liquefaction potential of a seabed foundation.

Earthquake-induced soil deformation is an important factor in the design of marine structures in the earthquake active regions. Numerous empirical or semi-empirical approaches have considered the influence of the geology, tectonic source, causative fault type, and frequency content of earthquake motion on lateral displacement caused by liquefaction. Pirhadi et al. [2] added an earthquake parameter of the standardized cumulative absolute velocity to the original dataset for analysis. They proposed a new response surface method (RSM) approach, which is applied on the basis of the artificial neural network (ANN) model to develop two new equations for the evaluation of the lateral displacement due to liquefaction.

#### **3. Foundations of Marine Infrastructures**

The stability of marine infrastructures is an important parameter in the design of offshore engineering projects. In this Special Issue, numerous marine infrastructures including pipelines, piled foundations, and submarine tunnels were investigated.

Offshore pipelines have been commonly used for the transportation of oil and gas. Therefore, safety of the pipeline route is one of the key factors in oil and gas projects. Unlike previous studies with FEM modeling, Wang et al. [6] proposed a meshfree model for the seabed, together with an OpenFOAM model for flow domain to examine the wave-induced transient soil response around an offshore pipeline. Both fully buried and partially buried pipelines in a trench layer were considered. Numerical examples demonstrated the capacity of their new meshfree model in the prediction of the wave-induced soil response. Foo et al. [7] adopted the FLOW-3D model together with poro-elastoplastic seabed model (within COMSOL) to examine the soil response around a fully buried pipeline under combined wave and current loading. They considered the residual soil response. In addition to wave and current loading, Zhang et al. [8] further considered earthquake loading for the wave–seabed–pipe interaction problem. In this study, they considered both oil pipe and gas pipe in the model and concluded that the difference between the two cases was minor.

Monopiles have been adopted as the supporting structures for various marine structures such as platforms, offshore wind turbine foundations, cross-sea bridge piers, etc. Liu et al. [9] conducted a series of laboratory tests for the dynamic response of offshore open-ended pile under lateral cyclic loading. They also used a discrete element model for numerical simulation and compared their results with the experimental data. Based on the numerical examples, they found that both the soil plug and outer friction contributed significantly to the pile lateral resistance; the "developing height" of the soil plug under lateral loading is in the range of two times the pile diameter above the pile end.

He et al. [11] employed ABAQUS to establish the interaction between rock-socketed monopile and layered soil–rock seabed. Based on a combined finite–infinite element model, the dynamic impedances and dynamic responses of large diameter rock-socketed monopiles under harmonic load are analyzed. When rock-socketed depth increases, the dynamic stiffness of pile increases, while the sensitivity to dimensionless frequency decreases. This indicates that the ability of pile to resist deformation increases under dynamic load, which is consistent with the results obtained from monopile deformation analysis.

In addition to geotechnical issues, the scouring of soil around large-diameter monopile will alter the stress history, and therefore the stiffness and strength of the soil at shallow depth, with important consequence to the lateral behavior of piles. The role of stress history was investigated for a larger diameter monopile [10]. Their study concluded that scour significantly increases the over-consolidation ratio and reduces the undrained shear strength of the remaining soil, which contributes to the significant difference in pile behavior between considering and ignoring the stress history effect.

Xiong et al. [12] developed a scour identification method, based on the ambient vibration measurements of superstructures. The Hangzhou Bay Bridge was selected to illustrate the application of the proposed model. Their study found that the high-order vibration modes are insensitive to the scour.

In addition to pipeline and pile foundation, based on COMSOL Multiphysics, Chen et al. [13] developed a two-dimensional coupling model of a wave–seabed–immersed tunnel for the dynamic responses of a trench under wave action in the immersing process of tunnel elements. Both liquefaction and shear failure are examined in this paper.

The buoyancy of the bottom-supported foundation is a critical issue in platform design because it counteracts parts of the vertical loads. In [14], a model box is designed and installed with earth pressure transducers and pore pressure transducers to simulate the sitting process of the bottom-supported foundation.

#### **4. Coastal Geohazard**

In this Special Issue, there are four papers related to other marine geohazard issues. Among these, Zhu et al. [15] reported the evidence of submarine landslide in South China Sea, and analyzed the causes of these events, based on their long-term field observations. Three concurrent events (the shoreward shift of the shelf break in the Baiyun Sag, the slump deposition, and the abrupt decrease in the accumulation rate on the lower continental slope) indicate that the giant Baiyun–Liwan submarine slide in the PRMB, South China Sea, occurred at 23–24 Ma, in the Oligocene–Miocene boundary. This landslide extends for over 250 km, with the total affected area of the slide up to 35,000–40,000 km2. Their research suggests that coeval events (the strike–slip movement along the Red River Fault and the ridge jump of the South China Sea) in the Oligocene–Miocene boundary triggered the Baiyun–Liwan submarine slide. Zhu et al. [16] developed a simple approach to investigate the stability of an unsaturated and multilayered coastal-embankment slope during the rainfall, in which a Random Search Algorithm (RSA) based on the random sampling idea of the Monte Carlo method was employed to obtain the most dangerous circular sliding surface, whereas the safety factor of the unsaturated slope was calculated by the modified Morgenstern–Price method. It was found that the fluctuation of the groundwater level has a significant influence on the location of the most dangerous sliding surface. The associated minimum safety factor and the sliding modes of unsaturated-soil slope gradually change from deep sliding to shallow sliding with the rise of groundwater level.

The stability of hydrate-bearing near-wellbore reservoirs is one of key issues in gas hydrate exploitation. A thermo-hydro-mechanical-chemical (THMC) multi-field coupling mathematical model considering damage of hydrate-bearing sediments is established in [17]. As reported in the paper, with continuous hydrate dissociation, the cementation of the sediment gradually decreases, and the structural damage gradually increases. This will lead to the partial softening and stress release of the stratum and will result in the decline of the bearing capacity of the reservoir. Therefore, damage of hydrate-bearing sediments has an adverse impact on the stability of the near-wellbore reservoir.

Li et al. [18] conducted a series of pumping well tests for the coastal micro-confined aquifer (MCA) in Shanghai to investigate the dewatering-induced groundwater fluctuations and stratum deformation. With the field tests, a numerical method is proposed for the estimation of hydraulic parameters and an empirical prediction method is developed for dewatering-induced ground settlement. The proposed prediction method worked well in most of the test site except in the far-field and the central parts. The parameters used in the method can be obtained by performing fitting with observation data, avoiding the dependence on precise hydrogeological parameters.

**Author Contributions:** D.-S.J., writing—original draft preparation; J.Z., writing—review and editing; and Ö.K., writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful for the supports from all authors and reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Wave (Current)-Induced Pore Pressure in Offshore Deposits: A Coupled Finite Element Model**

#### **Chencong Liao 1,\*, Dongsheng Jeng 1,2, Zaibin Lin 3, Yakun Guo <sup>4</sup> and Qi Zhang <sup>1</sup>**


Received: 10 May 2018; Accepted: 3 July 2018; Published: 6 July 2018

**Abstract:** The interaction between wave and offshore deposits is of great importance for the foundation design of marine installations. However, most previous investigations have been limited to connecting separated wave and seabed sub-models with an individual interface program that transfers loads from the wave model to the seabed model. This research presents a two-dimensional coupled approach to study both wave and seabed processes simultaneously in the same FEM (finite element method) program (COMSOL Multiphysics). In the present model, the progressive wave is generated using a momentum source maker combined with a steady current, while the seabed response is applied with the poro-elastoplastic theory. The information between the flow domain and soil deposits is strongly shared, leading to a comprehensive investigation of wave-seabed interaction. Several cases have been simulated to test the wave generation capability and to validate the soil model. The numerical results present fairly good predictions of wave generation and pore pressure within the seabed, indicating that the present coupled model is a sufficient numerical tool for estimation of wave-induced pore pressure.

**Keywords:** wave motion; offshore deposits; seabed response; FEM; pore pressure

#### **1. Introduction**

The phenomenon of wave and seabed interaction has drawn great interest among coastal geotechnical engineers over the past decade. The reason for this growing attention is that offshore infrastructure, such as platforms, pipelines and breakwaters, have encountered structural failure due to wave-induced seabed instability [1–3] rather than construction or material failure.

Considerable investigations into the wave-seabed interaction have been carried out in past decades. The methods for investigating the wave-seabed interaction problem mainly include three types, namely the uncoupled method, the semi-coupled method, and the fully coupled method [4–6]. The uncoupled method in investigating a wave-induced seabed response mainly occurred in earlier studies. There is no data exchange between the fluid motion and the seabed deformation. The porous and deformable seabed was regarded as a rigid and impermeable medium in a fluid domain, and the dynamic wave pressure on the seabed surface was replaced by a simplified wave pressure equation in the seabed domain [7,8]. The semi-coupled method, also called the one-way coupled method, has been widely used in investigating the wave-seabed interaction problem in past decades. The wave motion was firstly calculated through CFD (computational fluid dynamic) solver, which is usually coded by FDM (finite difference method) and FVM (finite volume method). Then, the dynamic response of the seabed was analyzed by FEM (finite element method), in which the dynamic wave pressure

extracted from the fluid domain was applied on the seabed surface through linear interpolation [9,10]. The semi-coupled method could consider the effects of the dynamic wave loading on the seabed. However, the feedback of the deformed seabed to the wave motion is not taken into account [11–14]. The fully coupled method could simulate both the wave motion and the dynamic seabed response simultaneously, in which a real-time data exchange is required between the two domains. It is easy to see that the fully coupled method should be the most accurate method for studying the wave-seabed interaction problem. However, in investigating the wave-seabed interaction problem, the fully coupled method is scarcely used in the previous research.

To implement the wave propagating process, it is necessary to build a wave-maker in the wave field, where the progressive wave is generated and propagates over a porous seabed combined with currents. Based on the FEM, we use an internal wave-maker method for generating essentially directional waves in a two-dimensional domain using a momentum source function of the Reynolds Averaged Navier Stokes (RANS) equation proposed by Choi and Yoon [15]. The internal wave-maker was used to avoid the influence of waves reflected from the wave-maker toward the domain because the waves generated by the source function do not interact with waves reflected from inside the domain and the sponge-layer method, as proposed by Israeli and Orszag [16], has been used to absorb outgoing waves generated by the wave-maker in the present study.

To sum up, both the wave and seabed field are modelled by FEM in this study. No interface program is needed to transfer the loads between them. The structure of the present paper is illustrated as follows. Section 2 introduces the basic equations that describe the wave-seabed interaction. The revised RANS equations govern the ocean wave, while the poro-elastoplastic equations describe the mechanical behavior of the seabed under wave loading. In Section 3, the present model is validated against the analytical solution and the available data of experiments shown in the literature. This section includes the wave module verification, seabed module verification and wave-seabed interaction verification. Finally, the application of the present model on wave-induced pore pressure and liquefaction is illustrated in Section 4.

#### **2. Theory and Methods**

Two sub-modules are included in the present coupled approach: A wave module and a seabed module, as shown in Figure 1. The wave module is established in order to generate the wave train (current) and to describe the viscous wave propagation. The seabed module is adopted to calculate the seabed response to wave loading. Unlike any previous one-way coupled models, both sub-modules are strongly integrated in COMSOL Multiphysics (COMSOL 5.2) [17].

**Figure 1.** Sketch of wave (current)-seabed interaction.

#### *2.1. Wave Module*

In the present study, the internal wave-maker [15] was adopted to generate a progressive wave with sponge layers [16] to absorb the wave at both ends of the numerical flume. Thus, the wave reflection from both flume ends could be efficiently eliminated.

The wave propagation above the porous seabed is described by solving the revised Reynolds Averaged Navier Stokes (RANS) equations, which are derived by integrating the momentum source term into the RANS equations, and which govern the wave motion in an incompressible fluid:

$$\frac{\partial u\_l}{\partial x\_i} = 0 \tag{1}$$

$$\frac{\partial u\_i}{\partial t} + u\_j \frac{\partial u\_i}{\partial u\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial x\_i} + g\_i + \frac{1}{\rho} \frac{\partial \tau\_{ij}}{\partial x\_j} \tag{2}$$

where *i, j* =1, 2, 3 denotes the dimensions of wave motion; *ui* is the *i*th component of fluid velocity; *ρ* is the fluid density; *p* is the fluid pressure; *gi* is the gravitational force; and *τij* is the viscous stress tensor.

The *k-ε* model is employed to enclose the turbulence:

$$\frac{\partial \mathbf{k}}{\partial t} + \mathbf{u}\_{\dot{f}} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\dot{f}}} = \frac{\partial}{\partial \mathbf{x}\_{\dot{f}}} \left[ \left( \frac{\nu\_{\text{f}}}{\sigma\_{\text{k}}} + \nu \right) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\dot{f}}} \right] + 2\nu\_{\text{f}} \sigma\_{\text{i}\dot{f}} \frac{\partial \mathbf{u}\_{\dot{i}}}{\partial \mathbf{x}\_{\dot{f}}} - \varepsilon \tag{3}$$

$$\frac{\partial \varepsilon}{\partial t} + u\_j \frac{\partial \varepsilon}{\partial x\_j} = \frac{\partial}{\partial x\_j} \left[ \left( \frac{\nu\_l}{\sigma\_\varepsilon} + \nu \right) \frac{\partial \varepsilon}{\partial x\_j} \right] + C\_{1\varepsilon} \frac{\varepsilon}{k} (2\nu\_l \sigma\_{i\bar{j}}) \frac{\partial u\_l}{\partial x\_{\bar{j}}} - C\_{2\varepsilon} \frac{\varepsilon^2}{k} \tag{4}$$

where *k* is the turbulent kinetic energy; *ν* is the kinematic viscosity; and *ν<sup>t</sup>* = *Cdk*2/*ε* is the eddy viscosity with *Cd* = 0.09. The empirical coefficients are *σ<sup>k</sup>* = 1.0, *σε* = 1.3, *C*1*<sup>ε</sup>* = 1.44 and *C*2*<sup>ε</sup>* = 1.92. The rate of stress tensor is *σij* = <sup>1</sup> 2 *∂ui ∂xj* + *∂uj ∂xi* , and *τij* = 2 *ν* + *Cd k*2 *ε <sup>σ</sup>ij* <sup>−</sup> <sup>2</sup> <sup>3</sup> *kδij*, where *δij* is the Kronecker delta.

Generally speaking, there are several options to numerically generate a target wave via an internal wave-maker: Adding a mass source term to the mass conservation equation (Equation (1)) or introducing a momentum one to the equation of momentum conservation (Equation (2)). One can also use both the mass and momentum source to generate a train of wave. Theoretically, this mass/momentum source could be a point, line, or a finite volume source [18]. In this study, we will only focus on the issue of generating waves taking the method of a momentum source function in a two-dimensional domain.

To generate a wave with a momentum source, Equations (1) and (2) should be revised as follows:

$$\frac{\partial u\_i}{\partial x\_i} = 0 \text{ in } \Omega \tag{5}$$

$$\frac{\partial u\_i}{\partial t} + u\_j \frac{\partial u\_i}{\partial x\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial x\_i} + g\_i + \frac{1}{\rho} \frac{\partial \mathbf{r}\_{ij}}{\partial x\_j} + S\_i \text{ in } \Omega \tag{6}$$

where *Si* is the momentum source function within a finite area Ω. Once the simulation starts, the free surface above the source region (Ω) will vibrate instantly and the surface vibration starts to propagate to both ends of the wave flume.

To properly explain the expression of *Si* in Equation (4), it is necessary to relate the mass source function to the momentum source function for wave generation:

$$S\_i = \left(\mathbf{S}\_{\mathbf{x}}, \mathbf{S}\_{\mathbf{y}}\right) = -\mathbf{g}\nabla \left(\int f \mathbf{d}t\right) \tag{7}$$

*J. Mar. Sci. Eng.* **2018**, *6*, 83

There are several expressions of *f* due to different wave generation theories with a mass source. The following expression is adopted from the revised Boussinesq's equation [19]:

$$f(\mathbf{x}, y, t) = \frac{\exp\left(-\beta \mathbf{x}^2\right)}{4\pi^2} \int\_{-\infty}^{+\infty} \int\_{-\infty}^{+\infty} D(k, \omega) \exp[i(ky - \omega t)] d\omega dk \tag{8}$$

where:

$$D = \frac{2A\_0(\omega^2 - a\_1gk^4h^3)\cos\theta}{\omega I\_1k[1 - a(kd)^2]}\tag{9}$$

in which the angular frequency, *ω*, water depth, *d*, wave number, *k*, wave obliquity, *θ*, and wave amplitude, *A*0, are the wave parameters adopted to obtain a target wave train. In addition, *<sup>I</sup>*<sup>1</sup> = *<sup>π</sup>*/*<sup>β</sup>* exp(−*k*2/4*β*), where *<sup>β</sup>* = 80/*δ*2/*L*2, in which *<sup>L</sup>* is the wavelength and *<sup>δ</sup>* is a parameter characterizing the width of the internal wave generation region. Another expression is from the revised RANS equations proposed by Lin and Liu [18] as follows:

$$\int\_{0}^{t} \int\_{\Omega} f(x, y, t) \mathrm{d}\Omega \mathrm{d}t = 2 \int\_{0}^{t} \mathrm{C}\eta(t) \mathrm{d}t \tag{10}$$

where *C* is the wave velocity and *η*(*t*) is the free surface elevation above the source region. By using adequate wave parameters in Equations (5), (6) and (8), any target wave can be obtained.

Regarding the process of current generation, the steady current flow is generated in the whole domain before wave generation. Once the current becomes stable, the internal wave maker starts to generate a wave. Then, the current and wave are coupled and the wave propagates from the wave-maker zone towards the sponge areas at both ends of water domain.

#### *2.2. Seabed Module*

The wave-induced pore pressure, *pe*, varies with time at a given location as suggested by Sassa and Sekiguchi [20], and consists of two components:

$$p\_{\varepsilon} = p\_{\varepsilon}^{(1)} + p\_{\varepsilon}^{(2)} \tag{11}$$

where *p* (1) *<sup>e</sup>* is oscillatory pore pressure and *p* (2) *<sup>e</sup>* represents the residual component.

#### 2.2.1. Oscillatory Response of Soil

On the basis of the conservation of mass equation, Biot's consolidation equation [21] are adopted as the governing equation for oscillatory response. For two-dimensional analysis, the mass conservation is expressed as follows:

$$\frac{\partial^2 p\_\varepsilon^{(1)} }{\partial x^2} + \frac{\partial^2 p\_\varepsilon^{(1)} }{\partial z^2} - \frac{\gamma\_w n\_s \beta\_s}{k\_s} \frac{\partial p\_\varepsilon^{(1)} }{\partial t} + \frac{\gamma\_w}{k\_s} \frac{\partial \varepsilon\_\varepsilon}{\partial t} = 0 \tag{12}$$

where *γ<sup>w</sup>* and *ns* denote the unit weight of water and the soil porosity, respectively. The volume strain, *εe*, and the compressibility of pore fluid, *βs*, are defined as, respectively:

$$
\varepsilon\_{\varepsilon} = \frac{\partial u\_s}{\partial x} + \frac{\partial w\_s}{\partial z} \text{ and } \beta\_s = \frac{1}{K\_{\text{uv}}} + \frac{1 - S}{P\_{\text{uv}}} \tag{13}
$$

where (*us*, *ws*) are the soil displacements; *Kw* is the true modulus of elasticity of pore water (taken as <sup>2</sup> <sup>×</sup> 109N/m2 [22]); *Pwo* is the absolute water pressure; and *<sup>S</sup>* is the degree of saturation.

The total stress, *σij*, can be decomposed into the effective stress, *σij* , and the pore pressure:

$$
\sigma\_{ij}^{(1)} = \sigma\_{ij}^{\prime(1)} + \delta\_{ij} p\_c^{(1)} \tag{14}
$$

where *δij* is the Kronecker delta.

Ignoring the body forces, the equilibrium equations can be written as:

$$
\delta G \nabla^2 u\_s + \frac{G}{(1 - 2\mu\_s)} \frac{\partial \varepsilon\_\varepsilon}{\partial x} = -\frac{\partial u\_c^{(1)}}{\partial x} \tag{15}
$$

$$
\delta G \nabla^2 w\_s + \frac{G}{(1 - 2\mu\_s)} \frac{\partial \varepsilon\_\varepsilon}{\partial z} = -\frac{\partial u\_\varepsilon^{(1)}}{\partial z} \tag{16}
$$

in the *x*− and *z*− directions, respectively.

Equations (10), (13), and (14) are the governing equations accounting for the oscillatory mechanism, in which the undetermined soil displacements and oscillatory pore pressure are to be solved.

In accordance with elastic theory, other stresses can be written, based on soil displacements, as:

$$\sigma\_x'^{(1)} = 2G \left[ \frac{\partial u\_\delta}{\partial x} + \frac{\mu\_\delta \varepsilon\_\varepsilon}{1 - 2\mu\_\delta} \right], \\ \sigma\_z'^{(1)} = 2G \left[ \frac{\partial w\_\delta}{\partial z} + \frac{\mu\_\delta \varepsilon\_\varepsilon}{1 - 2\mu\_\delta} \right], \\ \tau\_{xz}^{(1)} = G \left[ \frac{\partial u\_\delta}{\partial z} + \frac{\partial w\_\delta}{\partial x} \right] = \tau\_{zx}^{(1)} \tag{17}$$

#### 2.2.2. Residual Response of Soil

Following Sassa and Sekiguchi [20], Liao et al. [23] extended the one-dimensional model to a two-dimensional model. In the model, the evolution of the residual pore pressure, *p* (2) *<sup>e</sup>* , can be expressed as:

$$\frac{\partial p\_{\epsilon}^{(2)}}{\partial \xi^{\mathbb{F}}} = \frac{2\pi K\_{\upsilon}k\_{\ast}}{\omega \gamma\_{\upsilon\upsilon}} \left( \frac{\partial^{2} p\_{\epsilon}^{(2)}}{\partial x^{2}} + \frac{\partial^{2} p\_{\epsilon}^{(2)}}{\partial z^{2}} \right) + K\_{\upsilon} R \beta e^{-\beta \xi} (e^{a\chi} - 1) \tag{18}$$

where *Kv* = E/2(1 − *μs*) represents the bulk modulus of soil. The expression of plastic volumetric strain can be written as:

$$\varepsilon\_p(\xi,\chi) = \varepsilon\_p^{\infty}(\chi) \cdot \left[1 - e^{-\beta \tilde{\varepsilon}}\right], \varepsilon\_p^{\infty}(\chi) = R \cdot \left[e^{a\chi} - 1\right] \tag{19}$$

where *R* and *α*, *β* are the parameters of material. The cyclic stress ratio, *χ*, can be expressed as:

$$\chi(\mathbf{x}, z) = \frac{|\pi(\mathbf{x}, z)|}{\sigma\_{\mathbf{r}0}'(z)}\tag{20}$$

where *τ*(*x*, *z*) is the maximum amplitude of shear stress and *σ <sup>v</sup>*0(*z*) stands for the initial effective stress in the vertical direction.

In Equation (18), the first term on the right-hand side (RHS) represents the rate of pore pressure build-up and dissipation. The second term on the RHS correlates to the effect of cyclic loading (wave repetition) on the accumulation of residual pore pressure. For more details of the poro-elastoplastic model, the readers can refer to Sassa and Sekiguchi [20] and Liao et al. [23].

#### *2.3. Coupling Method*

In this section, the coupled process of the present model will be presented, including the time scheme, the mesh scheme, and the boundary conditions.

#### 2.3.1. Time Scheme

In the present study, the identical time scheme is applied to the whole computation domain. Since the seabed module easily reaches convergence, the time interval is set to be adaptive, thus fulfilling the requirement of fluid flow. In the traditional model, the non-matching time scheme may also work for the present case. However, it produces cumulative errors in interpolating time steps between the wave module and the seabed module. Furthermore, the non-matching time

scheme makes the information exchange between the two sub-modules more complicated. In the authors' opinion, the non-matching time scheme may reduce the accuracy of computation, thus the matching time scheme is applied to achieve a more accurate computation. FEM is used to solve all the governing equations, in combination with the second-order Lagrange elements, to ensure the second order of accuracy in evaluating the dependent variables in the computational domain. The Generalized-α Method was used for the time integration when computing the dynamic soil response under water action. As a second order accurate numerical scheme, the Generalized-α Method is a one-step, three-stage implicit method, in which accelerations, velocities and displacements are uncoupled. To obtain computational stability, the time interval is automatically adjusted to satisfy the Courant–Friedrichs–Lewy condition and the diffusive limit condition.

#### 2.3.2. Mesh Scheme

In the process of solving the revised RANS equations and elastoplastic equations, two typical mesh types (i.e., the matching mesh and the non-matching mesh) are adopted in the present study. The matching mesh requires the same numbers of mesh nodes along the seabed surface. However, the solid element size is generally much larger than that of the fluid cells to reach the acceptable computation efficiency. Therefore, it is necessary to use the non-matching mesh system outside the seabed surface to make sure each part of the models is calculated in proper meshes. This treatment of the mesh scheme will not affect the process of information exchange between the two modules because the matching mesh is applied at the seabed surface; particularly, the application of non-matching mesh helps reduce the cost of CPU time and memory occupation.

The meshes used in the fluid domain are structured four-node quadrilateral elements, and the simulated results are broadly affected by the grid resolution. As such, certain criteria should be satisfied to generate a high-quality mesh to ensure a valid, and hence accurate, solution. The model grid sensitivity studies show that the model is convergent using the resolution of mainly L/60 in the x- and y- directions and H/10 in the z- direction, with a refinement factor of 2, where L is the wave length; H is the wave height; and the refinement factor represents the ratio between the grid solution of the area without structural influence and the refinement area in the vicinity of the structure. The optimal non-orthogonal FEM meshes are used in the seabed domain, which are automatically generated by the COMSOL software with a maximum element size scaling factor (MESSF) controlling the maximum allowed element size. The mesh is refined until no significant changes in the numerical solution was achieved.

#### 2.3.3. Boundary Conditions

Appropriate boundary conditions are required to close the problem of wave-seabed interaction. Firstly, in the water domain the upper boundary of the air layer is set as a pressure outlet, where the pressure can flow in and out without any constraint. Secondly, the continuity of pressure and fluid displacement is applied at the air/water interface. Then, at the bottom boundary of the water domain, the displacement of the water particles is equal to that of the seabed surface.

Following the previous studies [24], it is acceptable to set the vertical effective stress and shear stresses to zero at the seabed surface:

$$\sigma\_z' = 0, \ \tau\_{xz} = \tau\_b(\mathbf{x}, t), \ p\_c^{(1)} = P\_b, \ p\_c^{(2)} = 0, \text{at} \ z = 0 \tag{21}$$

where *Pb*(*x*, *t*) and *τb*(*x*, *t*) are the wave pressure and shear stress at the seabed surface, respectively, and can be obtained from the wave model outlined in Section 2. Secondly, for the soil resting on an impermeable rigid bottom, zero displacements and no vertical flow occurs at the horizontal bottom:

$$u\_{\varepsilon} = w\_{\varepsilon} = 0, \frac{\partial p\_{\varepsilon}^{(1)}}{\partial z} = \frac{\partial p\_{\varepsilon}^{(2)}}{\partial z} = 0, \text{at } z = -h,\tag{22}$$

#### 2.3.4. Coupled Process

In the coupled process, the wave module is in charge of the simulation of the wave (current) propagation and determines the wave loading. The standard *k*-*e* turbulence model with the level set method (LSM) and the moving mesh method are used to model the flow of two different, immiscible fluids, where the exact position of the interface is of interest. The interface position is tracked by the LSM, with boundary conditions that account for surface tension and wetting, as well as mass transport across the interface. The LSM tracks the air-water interface using an auxiliary function. Since the displacement of the seabed surface from the seabed module will affect the flow field in the wave module, the authors use the moving mesh method to track the time-dependent displacement of the seabed surface as well.

The seabed is modeled with the PDE (partial differential equation) interface in COMSOL Multiphysics to solve all the equations describing the elastoplastic soil. The wave pressure and forces acting on the seabed are simulated by the wave module, and the results are sent to the seabed module to capture the seabed response, mainly the displacements, pore pressure, and the effective stresses. Meanwhile, the feedback of seabed response to the flow field is taken into account without any time lag, thus achieving the coupling effect. The seawater and seabed displacements at the water-seabed interface are set to be identical as well as the pressures of seawater and pore water in the seabed. This boundary condition is the basic and key requirement that ensures the coupling process stated in this research.

#### **3. Model Validation**

In this section, the coupled model is validated against the analytical solution and the available data of experiments in the literature. This section includes the wave module verification, seabed module verification, and wave-seabed interaction verification.

#### *3.1. Wave Verification: Comparison with an Analytical Solution*

To validate the wave module, the free surface elevation and the dynamic water pressure that acts on the seafloor from the coupled model is verified against an analytical solution in terms of free surface elevation and the wave pressure on the seabed. The analytical solution is calculated by the Airy wave theory and compared with our numerical solution. The input data are as follows: Wave period, *T* = 12.0 s; water depth, *d* = 30 m; wave length, *L* = 170 m; and wave amplitude *η* = 2.5 m.

To simplify the discussion, we only show the data from the left half of the area for the analysis because the area is symmetrical. As shown in Figure 2, the domain is separated into two parts: The wave-propagation area and the wave-absorbing area. The wave is generated from the inner wave-maker zone and propagates through the entire water domain, and then absorbed by the sponge layer settled in the absorbing region. Figure 2a shows the spatial distribution of the free surface elevation from the present model and the analytical solution. The results agree well with each other except that the wave generated is slightly higher in our model than in the analytical solution. This is because when the wave approaches the absorbing area, the minor reflection from the sponge layer may amplify the wave height. This phenomenon can be alleviated by a proper absorbing coefficient or by using the data from areas away from the sponge layer. Similar results can be found in Figure 2b. The dynamic water pressure that acts on the surface of the seabed is slightly higher in our model than in the analytical solution due to the wave reflection of the sponge layer. Overall, the proposed model agrees well with the analytical solution in both the free surface elevation and the water pressure on the seabed surface.

**Figure 2.** Comparison of (**a**) free surface elevation and (**b**) water pressure between the coupled model and the analytical solution. The small amplitude wave theory is applied for the analytical solution.

#### *3.2. Seabed Verification: Comparison with Experimental Data*

Both oscillatory and residual soil response are verified in this section. The oscillatory pore pressure is compared with a one-dimensional compressive test conducted by Liu et al. [25], while the residual pore pressure is compared with a centrifugal test under water waves.

#### 3.2.1. Validation of the Oscillatory Pore Pressure

Liu et al. [25] conducted a series of one-dimensional laboratory tests to explore the vertical distribution of pore pressure under wave loading. The cylinder consisted of 10 cylindrical organic glass cells, as shown in Figure 3. Ten pore pressure gauges were installed in the sandy deposit, while one more pressure gauge was installed at the surface of the seabed. As presented in their study, only the oscillatory mechanism of the pore pressure was observed. Thus, the authors compare the results of oscillatory pore pressure with the data from laboratory experiments [25]. The simulated results of the vertical distribution of the maximum oscillatory pore pressure (*p* (1) *<sup>e</sup>* /*pb*) are illustrated in Figure 4.

**Figure 3.** Schematic diagram of one-dimensional cylinder equipment (adapted from [25]).

**Figure 4.** Comparison of oscillatory pore pressure with one-dimensional experimental data. The experimental data include pore pressure records of ten gauges (P0 to P9).

It should be noted that in the laboratory experiment, only the one-dimensional cylinder model facility was used. Thus, the wave length should be revised as infinite in the present model. Other input data used are included in Figure 4. The present model reaches a good agreement with the data from the one-dimensional experiment in the upper zone of the sandy soil. However, some discrepancy is observed in the deeper zone. There are two possible reasons that account for this discrepancy. The first is that the thickness of the sediments varies with the wave loading, which induces changes in the relative depth of sandy soil, resulting in pore pressure differences at the maximum amplitudes. Since the formulation used in the numerical approximation is based on elastic theory, accurate evaluation of the dynamic process with large soil deformation is not possible. Another potential reason is the variation of the soil density with depth in the experiment, since the sandy deposit is thick. Therefore, the response of the soil from the experimental tests may differ from that of the numerical curve, which was derived by assuming that soil properties are constant along the soil depth. However, the good agreement between the coupled model and experimental results is promising for prediction of the oscillatory pore pressure by the coupled approach.

#### 3.2.2. Validation of the Residual Pore Pressure

Sekiguchi et al. [26] conducted the first centrifugal standing wave tests to study the instability of horizontal sand deposits by a centrifugal method. A cross section through the wave tank is shown in Figure 5. The wave tank consisted of a wave channel, a wave paddle, and a sediment trench. Standing waves were formed under the condition of a frequency, *f* = 8.8 Hz, under steady 50 g acceleration, along with a fluid depth of 47 mm. Waves corresponded to a prototype condition of *d* = 2.35 m and *f* = 0.176 Hz. The amplitude of the input pressure,*p*0, was 1.7 kPa. The plastic parameter, *β*, was taken as 1.4 (corresponding to the parameter *α* in Sekiguchi et al. [26]). The parameters were set as follows: *<sup>α</sup>* = 55, *<sup>R</sup>* = 1.8 × <sup>10</sup><sup>−</sup>5, and *ns* = 0.5. More detailed information can be found in Sekiguchi et al. [26]. The excess pore pressure response measured in the centrifugal test is now compared with the prediction from the present poro-elastoplastic solution. As illustrated in Figure 6, the maximum pore pressure of the centrifugal test is slightly larger than that predicted by our model. This may be due to the fact that the water waves used in the present model would attenuate over the porous seabed, leading to deviation of the water pressure from the genuine value. Except for this, the simulation reaches a good agreement with the centrifugal test.

**Figure 5.** Cross-section of two-dimensional centrifuge equipment (units are mm).

**Figure 6.** Comparison of residual pore pressure with standing wave centrifugal test data [20]. Time history of pore pressure is from the soil element at elevation z/h = 0.25.

#### *3.3. Wave-Seabed Interaction Verification*

To the best of our knowledge, no laboratory experiments have been carried out to explore the interaction between the wave (or current) with the seabed except for Qi and Gao [27]. In their experiment, a series of laboratory tests was performed within a flume for the scour development and pore pressure response around a mono-pile foundation under the effect of combined wave and current. Wave, current, seabed, and structure were integrated in one model to investigate the interaction process, which provided a comprehensive understanding of the coupled model. In their experiments, the flow velocities in the wave field and pore pressure response around the pile in a finite seabed were measured simultaneously (Figure 7). Herein, the measured flow velocity and pore pressure of points far away from the mono-pile foundation, which represent the case of the wave (current)-seabed interaction without a structure, are selected for comparison with the present solution. The wave and current parameters in their experiment were: Water depth, *d* = 0.5 m; wave period, *T* = 1.4 s; wave height, *H* = 0.12 m; and current velocity, *U0* = −0.1, 0, and +0.1 m/s. The properties of the soil provided in their paper were: Shear modulus, *<sup>G</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>7N/m2 ; Poisson's ratio,*u* = 0.3; permeability, *<sup>K</sup>* = 1.88 × <sup>10</sup>4m/s; the void ratio, *<sup>e</sup>* = 0.771; and the soil was almost fully saturated.

**Figure 7.** Schematic diagram of three-dimensional water flume system (adapted from [27]). (Units are m).

The results with various current velocities are shown in Figures 8–10. The flow velocity represents the fluid particle velocity located at 0.2 m above the seabed surface. Differences between the present model and the experimental data are observed because the wave height generated in the experiment cannot be exactly the same as the value expected. Furthermore, the discrepancy occurs close to the wave crest and trough, which may be induced by the transformation of the linear wave profile into a nonlinear wave during propagation from the wave generator to the flume end. The phase of the pore pressure in porous sediments closely corresponds to the phase of the progressive wave above it. In spite of this, there is a trend toward overall agreement with the experimental data. Note, that in these cases, both the water wave (with current) and the seabed response are fully integrated into coupled FEM codes to simulate the wave (current)-seabed interaction, demonstrating the efficiency and ability of the present model when estimating the pore pressure in complex and multi-phase marine deposits.

**Figure 8.** Comparison of (**a**) flow velocity and (**b**) pore pressure between the present coupled model and experimental data (U0 = −0.1 m/s). The flow velocity represents the fluid velocity at 0.2 m above the seabed surface and pore pressure is from the soil element that was located at point C in Figure 7.

**Figure 9.** Comparison of (**a**) flow velocity and (**b**) pore pressure between the present coupled model and experimental data (U0 = 0 m/s). The flow velocity represents the fluid velocity at 0.2 m above the seabed surface and pore pressure is from the soil element that was located at point C in Figure 7.

**Figure 10.** Comparison of (**a**) flow velocity and (**b**) pore pressure between the present coupled model and experimental data (U0 = +0.1 m/s). The flow velocity represents the fluid velocity at 0.2 m above the seabed surface and pore pressure is from the soil element that was located at point C in Figure 7.

#### **4. Model Application**

Under cyclic wave loading, pore pressure varies extensively in the seabed. When wave-induced pore pressure exceeds a certain limit, soil liquefaction occurs. The following liquefaction criterion could be used to estimate the liquefaction potential:

$$\frac{1}{3}(\gamma\_s - \gamma\_w)(1 + 2K\_0)z \le p\_\mathfrak{c}(\mathbf{x}, y, z) - p\_\mathfrak{b}(\mathbf{x}, y) \tag{23}$$

where *γ<sup>s</sup>* and *γ<sup>w</sup>* are the unit weights of the seabed soil and water, respectively, and *pe*(*x*, *y*, *z*) and *pb*(*x*, *y*) are the wave-induced pore pressure in the seabed and the wave-induced pressure on the seabed surface, respectively.

To examine the difference between the coupled and uncoupled models on the wave-induced soil response in marine sediments, the development of wave-induced pore pressure and liquefaction depth is presented in Figure 11. To control other variables (like plasticity) that may affect the liquefaction area, both models only apply the elastic theory to modelling the seabed response. Other parameters are shown in Table 1.

**Figure 11.** Comparison of (**a**) pore pressure and (**b**) liquefaction depth within the seabed between coupled and uncoupled models.


**Table 1.** Input data for application of present model.

As shown in Figure 11, both wave-induced pore pressure and liquefaction depth in the coupled model are smaller than that in the uncoupled or semi-coupled model. Considering that the seabed surface displacement induced by water waves may alter the flow field in the wave generation model, which was not considered in the previous one-way model, the water pressure acting on the seabed suffers decreases due to the seabed surface motion. Therefore, the soil response may be slightly smaller than in the uncoupled model, and so may the liquefaction area. The results imply that the previous one-way or semi-coupled models ignored the attenuation effect of the porous seabed on water waves, resulting in an over-estimation of the wave-induced pore pressure within marine sediments. Otherwise, a physical scale model would be necessary to verify the results from the traditional uncoupled numerical model. Although, in this case, the discrepancy did not seem significant, its effects may be amplified when evaluating seabed liquefaction in the vicinity of marine structures. This conclusion would be of practical value when applying the traditional one-way model to evaluate the soil response during water wave loading.

#### **5. Conclusions**

In this paper, a coupled research method for solving wave-seabed interaction problem was presented. Revised RANS equations were employed to govern the ocean wave and the porous fluid in the seabed, while poro-elastoplastic equations were used to describe the mechanical response of the seabed under dynamic wave loading. The coupled numerical model was validated by comparison with an analytical solution (wave module) and experimental data (seabed module) in the literature. Overall, the consistency between the proposed model and the experimental data illustrated the capacity of predicting the response of the soil due to wave loading.

The major advantages of the coupled model include: (1) The wave and seabed models are coupled in the same platform (COMSOL Multiphysics) and all the equations are solved simultaneously; (2) the wave model can be used to simulate not only small amplitude waves but also large waves and nonlinear waves; (3) the elastoplastic model may be more precise when the plasticity of soil cannot be ignored, which is usually important for offshore deposits; and (4) the coupled model could be used for more complex situations such as the wave-seabed-structure interaction.

This paper presented the basic theory of a coupling model and compared it with the data available in the literature. The wave-seabed interaction with a marine structure, such as a pipeline or breakwater, could also be simulated with the present model. Research that focuses on the wave-seabed-structure issue will be carried out in the future.

**Author Contributions:** C.L. conducted the numerical analysis and write the paper. D.J. and Y.G. motivated this study with good and workable idea. Z.L. helped to build the wave maker in the software. And Q.Z. did the literature review for this paper.

**Acknowledgments:** The authors are grateful for the support from National Natural Science Foundation of China (Grant Nos. 41602282, 51678360 and 41727802)

**Conflicts of Interest:** The authors declare no conflicts of interest

#### **Nomenclature**



#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **New Equations to Evaluate Lateral Displacement Caused by Liquefaction Using the Response Surface Method**

#### **Nima Pirhadi, Xiaowei Tang and Qing Yang \***

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China; nima.pirhadi@yahoo.com (N.P.); tangxw@dlut.edu.cn (X.T.)

**\*** Correspondence: qyang@dlut.edu.cn; Tel.: +86-411-84707609

Received: 8 January 2019; Accepted: 31 January 2019; Published: 4 February 2019

**Abstract:** Few empirical and semi-empirical approaches have considered the influence of the geology, tectonic source, causative fault type, and frequency content of earthquake motion on lateral displacement caused by liquefaction (*DH*). This paper aims to address this gap in the literature by adding an earthquake parameter of the standardized cumulative absolute velocity (CAV5) to the original dataset for analyzing. Furthermore, the complex influence of fine content in the liquefiable layer (*F*15) is analyzed by deriving two different equations: the first one is for the whole range of parameters, and the second one is for a limited range of *F*<sup>15</sup> values under 28% in order to the *F*15's critical value presented in literature. The new response surface method (RSM) approach is applied on the basis of the artificial neural network (ANN) model to develop two new equations. Moreover, to illustrate the capability and efficiency of the developed models, the results of the RSM models are examined by comparing them with an additional three available models using data from the Chi-Chi earthquake sites that were not used for developing the models in this study. In conclusion, the RSM provides a capable tool to evaluate the liquefaction phenomenon, and the results fully justify the complex effect of different values of *F*15.

**Keywords:** liquefaction; lateral displacement; response surface method (RSM); artificial neural network (ANN)

#### **1. Introduction**

When, during earthquake motion, pore water pressure rises because of applied dynamic loads, the loose saturated sand layer that is relatively close to the ground surface is liquefied. Liquefaction can be discovered through manifestations such as (1) a sharp decrease in the frequency content of a sand layer, (2) settlement, (3) flow slides, (4) sand boiling, (5) foundation failure, and (6) lateral displacement (*DH*).

The movements of sand blocks, which have destroyed and affected constructions and infrastructure, ranging from a few centimeters to some meters [1], have been reported. Lateral displacement can be significantly damaging for piles, piers, and pipe lines during and for a short time after earthquakes and causes more damage to structures and infrastructures than any other type of liquefaction-induced ground failure. In this phenomenon, the large blocks of soil move towards the free face or along the slope. Researchers have developed several different models and approaches to predict the *DH* caused by liquefaction for some decades. Some of them have proposed numerical approaches [2–6] such as the finite element method (FEM) and the finite difference method (FDM). Next to that, analytical approaches have been developed, for example, minimum potential energy [7] and the sliding block model [8–12].

Among them, due to the complicated input model parameters and difficulties in their calculations, as well as because of the complex mechanism of liquefaction, empirical and semi-empirical models are the most common models that have been performed and developed by engineers and researchers [13–17]. However, in most cases, because of the scarcity and shortage in their database, some aspects of this phenomenon, such as geology, fault type, and the effect of near-fault sites, have been ignored, with the exception of Zhang et al. [18], who used Japanese spectral attenuation models, or Bardet et al. [19], who considered peak ground velocity (*PGV*) to overcome this shortage and improved the model proposed by Youd et al. [20]. Nevertheless, a shortage of studies in geology and motion frequency effects still exists.

In 2006, Kramer [21] reported the result of substantial research on around 300 ground motion parameters and declared that the most efficient and sufficient intensity measure on liquefaction is one standardized form of cumulative absolute velocity (*CAV*), which eliminated amplitudes less than 5 cm/sec2 and is defined as *CAV*5. Sufficiency defines which parameter is independent to estimate the target (increasing pore water pressure herein), and efficiency expresses which parameter is able to predict the target with lower uncertainties [22]. This parameter quantifies aspects of applied frequency load, which can be affected by the near-fault region aspect and causative fault type of earthquakes. Hui et al. [23] proposed an index of *PGV* to peak ground acceleration (*PGA*) to characterize the effect of liquefaction on the piles in near-fault zones. Further, Kwang et al. [24], through performing some uniform cyclic simple shear laboratory tests, demonstrated that *CAV*<sup>5</sup> provides the highest correlation with *DH* among ground motion parameters. While the significant correlation between *CAV*<sup>5</sup> and the evaluation of liquefaction have been characterized, no attempt has been made to take it into the account when developing empirical and semi-empirical models.

Furthermore, artificial intelligence has been applied to develop models and correlations to predict *DH* using databases that were collected from sites [25–28]. Training is organized to minimize the mean square error (MSE) function. Wang et al. [25] used a back-propagation neural network to develop a model for the prediction of lateral ground displacements caused by liquefaction. They applied the same records used by Bartlet et al. [1], along with 19 datasets of Ambraseys et al. [29]. Among all datasets, 367 data points were used for the training phase, and the extra 99 datasets were used for the testing phase, while no validating phase was conducted. The model was developed using the same parameters suggested by Youd et al. [14].

Baziar et al. [26] created two subsets for training, to train a network to predict *DH*, and a validating phase, to prevent overtraining of the artificial neural network (ANN) model. Then, they presented an ANN model using STATISTICA software (version of Statistica 5.1, Dell Software, Round Rock, TX, USA) to estimate *DH*. They inspected the performance of their model using validating subset data without considering extra available models. Furthermore, a new model was presented by Javadi et al. on the basis of genetic programming (GP) [27]. They divided the dataset randomly, without paying attention to the statistical properties of the input parameters, into two subsets for the validating and training phase. Garcia et al. established a neuro-fuzzy model to use the advantages of both systems. They randomly separated their dataset into two subsets for training and testing; however, they did not take statistical aspects into account. They also compared the value predicted by their model with extra models to evaluate its performance. Baziar et al. [28] then applied ANN and GP to propose a new model. They divided their dataset randomly into two subsets for the testing and training phase; a validating process was not performed, and the statistical factors of the parameters were not considered.

Although the effects of fine content (*Fc*) in different values on excess pore water pressure have been investigated [30–34], to the best knowledge of the authors, no attempts have been made to consider the range of *Fc* to establish the models to predict *DH*. Most of the studies reveal a range of 20% to 30% for the transition of behavior of the response of sand to earthquake and liquefaction occurrence. Maurer et al. [35] investigated the Canterbury earthquakes in 2010 and 2011 through 7000 case history datasets and illustrated that a high value of *Fc* caused more inaccuracy in liquefaction assessments. Tao performed some laboratory tests and demonstrated that the potential of liquefaction has a significant dependency on initial relative density (*Dr*) when the *Fc* value is larger than 28% [32].

This study is based on the database of Youd et al. [20] and the addition of a new earthquake parameter of *CAV*5, which is *CAV* with a 5-cm/sec2 threshold acceleration, through the attenuation equation presented by Kramer et al. [21]. By adding *CAV*5, the dataset was expanded and became more capable of and efficient in considering aspects of earthquakes and geology site situations, such as earthquake motion frequency, near-fault effects and the causative fault type of an earthquake. The second dataset was created by eliminating samples with an average *Fc* in a liquefiable soil layer (*F*15) less than 28%. The response surface method (RSM) is used for the first time as a novel method to develop two equations to predict lateral displacement due to liquefaction (*DH*) in order to two created datasets herein. Furthermore, the meaningful and effective terms of the equations are discovered through hypothesis testing of the *p*-value. In this study, two ANNs with back-propagation analysis were developed to measure the coding input data of the RSM. To develop each ANN model, the main dataset is first divided into three subsets for the training, testing, and validating stages, considering statistical properties—instead of random division—to increase the capability and accuracy of the model. To achieve this goal, an attempt is made to create all three subsets with close statistical factors. Finally, the results are compared with data measured from the Chi-Chi earthquake's near fault zone of Wufeng district (Figure 1) and Nantou district (Figure 2), as well as with the predicted *DH* through three extra models [20,27,36] to demonstrate the accuracy and capability of the RSM models.

**Figure 1.** Location of lateral displacement and ground failure due to liquefaction during the Chi-Chi earthquake in the Wufeng district.

**Figure 2.** Location of lateral displacement and ground failure due to liquefaction during the Chi-Chi earthquake in the Nantou district.

#### **2. Review of Empirical and Semi-Empirical Models**

Bartlett and Youd [1], based on factors in References [13,14,29], developed a new model to predict *DH* due to liquefaction; they supposed that earthquake, topographical, geological, and soil factors are the most influential parameters on *DH*. They studied 467 displacement vectors from the case history database. Among those vectors, 337 were from the 1964 Niigata and 1983 Nihonkai-Chubu, Japan, earthquakes; 111 were from earthquakes in the United States; and the other 19 cases were selected from Ambraseys' [29] database. In the end, they developed a new model by using multiple linear regression (MLR) for free-face and ground slope conditions [37], but they did not separate earthquakes according to their region because of a database shortage. Youd et al. revised their MLR by adding case history data from three earthquakes (1983 Borah Peak, Idaho; 1989 Loma Prieta; and 1995 Hyogoken-Nanbu (Kobe)), and they considered coarser-grained materials. They removed eight displacement sites with prevented free lateral movement and developed two equations with more accuracy given as follows:

Free-face conditions:

$$\begin{array}{l} \log D\_{H} = & -16.713 + 1.532M \\ & -1.4\log r^\* - 0.012r + 0.592\log W + 0.540\log T\_{15} + 3.413\log(100 - F\_{15}) \\ & -0.795\log(D50\_{15} + 0.1\text{ mm}) \end{array} \tag{1}$$

Sloping ground conditions:

$$\begin{aligned} \log D\_H &= -16.213 + 1.532 M\_w - 1.406 \log r^\* - 0.012R + 0.338 \log S + 0.540 \log T\_{15} \\ &+ 3.413 \log(100 - F\_{15}) - 0.795 \log (D50\_{15} + 0.1 \text{ mm}) \end{aligned} \tag{2}$$

where

$$r^\* = r + r\_0 \tag{3}$$

and

$$r\_0 = 10^{(0.89M - 5.64)} \tag{4}$$

In Equations (1) and (2), *DH* is the predicted lateral ground displacement (m), *Mw* is the moment magnitude of the earthquake, and *T*<sup>15</sup> is the cumulative thickness of saturated granular layers (m) with corrected blow counts ((*N*1)60) less than 15. Moreover, *F*<sup>15</sup> is the average fines content of sediment within *T*<sup>15</sup> (%); *D*5015 is the average mean grain size for granular materials within *T*<sup>15</sup> (mm); *S* is the ground slope (%); and *W* is the free-face ratio (*H*/*L*), where *H* is the height of the free face and *L* is the distance from the base of the free face to the liquefied point. Finally, *r* is the nearest horizontal or map distance from the site to the seismic energy source (Km).

Rezania et al. [37] developed a model, based on evolutionary polynomial regression, for the assessment of liquefaction potential and lateral spreading. According to response spectral acceleration, measured through strong-motion attenuation models, Zhang et al. [38] revised the empirical model of Youd et al. [20] and demonstrated the ability of their model by comparing the predicted results with datasets from sites in Turkey and New Zealand [18]. Goh et al. proposed multivariate adaptive regression splines (MARS) by using data of Youd et al. [20]. They demonstrated an improvement of the original model [39].

#### **3. Artificial Neural Network**

An ANN is a computational process using a biological neural network structure. McCulloch et al. [40] were the first to introduce some simulations according to their neurology knowledge. Neural networks are strong approximators due to their ability to learn by samples and their independency from any algorithm or knowledge about internal features of the issue. Artificial neural networks have a number of advantages such as high accuracy in nonlinear relationships or dynamic mechanisms based on the number and effectiveness of samples. Neural networks are classically constructed in three types of layers. Layers are provided by interconnected nodes. Outlines of the network are developed via the input layer, which communicates to one or more hidden layers by performing a weighting process. Then, hidden layers link to a target (output layer). Further, learning is a supervised process that occurs with each cycle or "epoch" (i.e., each time the network is presented with a new input pattern).

Rumelhart et al. [41] introduced a back-propagation algorithm to decrease error according to the training data. The training process started with random weights to achieve minimum error. The calculation of the derivatives flows backwards through the network, which is why it is called back-propagation. The most common measure of error is the MSE:

$$\text{MSE} = \text{Ave } \| \text{(actual output vector } - \text{desired output vector)}^2 \| \tag{5}$$

Overtraining of a neural network happens when the network trains exactly to reply to just one kind of input, which is similar to rote memorization. Therefore, learning does not occur anymore. The ANN is able to be used for problems with non-linear or dynamic correlation.

#### **4. Response Surface Method**

The RSM is a group of statistical and mathematical techniques to develop a capable function for a relationship of response (*y*) or output variable, and input variables (*x*), given by:

$$y = f'(\mathfrak{x})\beta + \mathfrak{s} \tag{6}$$

where *f(x)* is a vector function, consisting of powers and cross-products of input variables. This function depends on the supposed form of the response.

There are some common forms that have been used by researchers. A second-degree polynomial with cross terms is the most complicated and the strictest among them, and it is given by:

$$R\_{(X)} = a\_0 + \sum\_{i=1}^n b\_i X\_i + \sum\_{i=1}^n \mathbb{C}\_i X\_i^2 + \sum\_{\substack{i=1 \\ i \neq j}}^n d\_{ij} X\_i X\_j \tag{7}$$

The main applications of RSM are defined as follows:


#### *4.1. Design of Experiments*

When more than one input factor is suspected of influencing an output, in order to fit physical or numerical experiments, a process by the name of design of experiments (DOE) [42] is developed. The DOE involves selecting some points according to which a response should be calculated.

In this paper, the design introduced by Box et al. [43] is used. It requires three levels to run an experiment. Furthermore, it is a special three-level design without any points at the vertices of the experiment region. This could be advantageous when the points on the corners of the cube represent level combinations that are prohibitively expensive or impossible to test because of physical process constraints. The design is applied in this study to prevent the input parameters' values from being negative. This DOE requires three levels of each input variable −1, 0, 1 (coded values) corresponding to minimum, middle, and maximum values of input parameters, respectively.

#### *4.2. Hypothesis Test*

The process in statistics science to meaningfully examine results is called a hypothesis test. During hypothesis tests, the validity of a claim, which is constructed about a population, is evaluated. This claim that is in essence on trial is called the null hypothesis. Hypothesis testing can be expressed in three steps:


One of the main approaches to make a decision to accept or reject a null hypothesis is the *p*-value, defined through an *α* value (the probability of error is called alpha). If the *p*-value is smaller than *α*, then the null hypothesis is rejected, and contrarily, if it is larger, then the null hypothesis is accepted.

In this paper, the common value of 0.05, which researchers have used for *α*, is applied to analyze the meaningfulness and significance of parameters and the terms in the RSM response (equation).

#### **5. Model Proposed**

Figure 3 illustrates the flow chart of the applied approach to develop the RSM equations to estimate *DH* in this study. Two models are presented: first one considered the whole range of the parameters and the second one was on the basis of the *F*<sup>15</sup> value being less than 18%.

**Figure 3.** Flowchart of the approach applied in this study to present an equation to predict lateral displacement caused by liquefaction.

#### *5.1. Dataset*

Bartlett et al. [37] collected 467 data samples of *DH* from the following eight earthquakes in the United States and Japan: San Francisco 1906, Alaska 1964, Niigata 1964, San Fernando 1971, Imperial Valley 1979, Borah Peak 1983, Nihonkai-Chubu 1983, and Superstition Hills 1987. The parameters of the case histories that they collected to analyze were divided into three groups:


Later, Youd et al. [20] eliminated eight sites' data from their dataset due to a lack of free lateral movement. Additionally, they added data from the following three earthquake sites: Borah Peak 1983, Loma Prieta 1989, and Hyogo-Ken Nanbu (Kobe) 1995.

In the present study, the main dataset was developed by adding a new parameter of *CAV*5, which is defined in Section 5.1.1. Kramer et al. [21] stated that *CAV*<sup>5</sup> is the most efficient and sufficient earthquake intensity to evaluate liquefaction in sandy soil. Therefore, *CAV*<sup>5</sup> was estimated using the attenuation equation presented by them.. In this way, the causative earthquake fault types of all earthquakes in the dataset were discovered. Then, the sites with a moment magnitude range from 6.4 to 7.9 were selected in order to find the applicable magnitude range for Equation (11); the Alaska 1964 site, with a magnitude of 9.2, was thus deleted from the dataset.

Furthermore, a statistical analysis was performed to examine and estimate the coefficient of correlation (*R*) of all input parameters with the output (*DH*) one by one. The estimated values of *R* for *r-DH* and *S-DH* were a positive value of 0.104 and a negative value of −0.98, respectively, contrary to the supposition for them. This was possibly due to the scarcity of sites that were explored, and

consequently, the shortage measured values for *r* and *S* in the main dataset. Therefore, in this study, after eliminating *r* and *S* from dataset, only the free-face condition was considered, and samples of the ground slope condition were deleted from the main dataset. Furthermore, *CAV*<sup>5</sup> was added to the dataset instead of *r*. Figures 4 and 5 plot *r* versus *CAV*<sup>5</sup> for range of *Mw* in the main dataset from 6.4 to 7 and 7 to 7.9, respectively. It should be mentioned that the bold points show more than one point coincided together.

**Figure 4.** *r* versus *CAV*<sup>5</sup> for 6.4 ≤ *Mw* < 7 of the main dataset applied in this study.

**Figure 5.** *r* versus *CAV*<sup>5</sup> for 7 ≤ *Mw* ≤ 7.9 of the main dataset applied in this study.

A dataset including 215 case histories with six parameters was then prepared. In addition, to investigate the complicated influence of fine content on the liquefaction, the second dataset was arranged by eliminating samples with an *F*<sup>15</sup> value larger than 28%. Therefore, the second dataset included 182 samples.

#### 5.1.1. Cumulative Absolute Velocity

Eed et al. utilized *CAV* as a criterion to evaluate the onset of structural damage for the first time. They reported it in the Electric Power Research Institute journal [44] and defined it in a mathematical framework as presented below:

$$CAV = \int\_0^{t\_{\text{max}}} |a(t)| dt \tag{8}$$

where *a*(*t*) is the acceleration of ground motion graph, *t* is the time, and *tmax* is the duration of the earthquake.

Liyanapathirana et al. [45] studied special aspects of Australian earthquakes and found that the predominant frequency of earthquakes in Australia is much higher than in California. This is because the earthquakes in Australia are in the middle of tectonic plates, so liquefaction has not occurred. They introduced a pseudo-velocity in the same dimensions and form as *CAV* to quantify this difference as indicated below:

$$V = \int\_0^{t\_{\max}} |a(t)| dt \tag{9}$$

By studying various Japanese codes, Orense demonstrated that seismic-induced shear stresses, calculated using empirical models containing *PGA*, do not have high accuracy in near-source circumstances but still have reasonable accuracy for far-source earthquakes. Next to that, he indicated threshold values of 150 gal and 20 kine for *PGA* and *PGV* respectively for the occurrence of liquefaction [46].

After inspecting approximately 300 parameters of earthquakes, Kramer et al. revealed that *CAV*5, as can be estimated through Equation (11), has better efficiency and sufficiency than other earthquake intensity parameters for liquefaction evaluation. They also utilized the strong Pacific Earthquake Engineering Research (PEER) database, consisting of 282 ground motions from 40 earthquakes to present an equation to calculate *CAV*<sup>5</sup> for shallow crustal events [21].

$$
\triangle A V\_{\sf S} = \int\_0^\infty \langle \mathfrak{X} \rangle |a(t)| \,\mathrm{where } \langle \mathfrak{X} \rangle = \begin{cases} \ 0 \ for \, |a(t)| < 5 \,\mathrm{cm}/\,\mathrm{sec}^2 \\\ \ 1 \ for \, |a(t)| \ge 5 \,\mathrm{cm}/\,\mathrm{sec}^2 \end{cases} \tag{10}
$$

$$
\ln\text{CAV}\_5 = 3.495 + 2.764(M - 6) + 8.539\ln(M/6) + 1.008\ln\left(\sqrt{r^2 + 6.155}\right) + 0.464\_1F\_N + 0.165F\_R \tag{11}
$$

where *CAV*<sup>5</sup> is a form of CAV based on Equation (10) (m/sec), *M* is the moment magnitude, and *r* is the closest distance to the rupture (km). *FN* = *FR* = 0 for strike slip faults, *FN* = 1 and *FR* = 0 for normal faults, and *FN* = 0 and *FR* = 1 for reverse or reverse-oblique faults.

The database of Equation (11) has a range of 4.7 to 7.4 for *M*, which is proposed for use for a maximum magnitude of 8, and it has a range of 1 to more than 100 km for *r*. In the present study, Equation (11) was applied to calculate *CAV*<sup>5</sup> from the dataset by considering the causative fault type.

#### *5.2. Artificial Neural Network Models*

At first, both datasets were divided into three groups: approximately 70% for the training phase and two groups of 15% each for the testing and validating phase. The validating phase was applied to prevent the model from being overtrained. The data deviation was conducted according to statistical factors, and the three groups contained similar statistical factors, such as minimum, maximum, and mean values. Furthermore, to achieve a higher accuracy of output, the same portion of any earthquake's data was selected for the three phases, so all the sites contributed with their data in the training, testing, and validating phases with similar statistical factors. The ANN was established using a back-propagation algorithm with one hidden layer, which is the most commonly used network to drive an ANN model to predict *DH*. There are six inputs, including six parameters, which were considered by Youd et al. [20] and Bartlett et al. [37]. In addition, there is the new measured parameter of *CAV*<sup>5</sup> as well as one output as a *DH*.

The first dataset was divided into three groups including 151 samples for training, 32 samples for testing, and an extra 32 samples for the validating phase. Tables 1 and 2 summarize the characteristics and certificates of the first dataset and the subsequently developed ANN model. As can be seen from Table 2, the coefficient of correlation (*R*) given in Equation (12), which is the most common factor to assess the performance of correlations, of the first ANN model for all three groups and the main dataset was around 90%.

$$R = \frac{\sum\_{i=1}^{N} (x\_i - x\_0)(d\_i - d\_0)}{\sqrt{\sum\_{i=1}^{N} \left(x\_i - x\_0\right)^2 \sum\_{i=1}^{N} \left(d\_i - d\_0\right)^2}}\tag{12}$$

where *n* is the sample size, *xi* and *di* are the individual sample points indexed with *i*, and *x* and *d* are the mean sample sizes.


**Table 1.** Characteristics of whole case histories' input parameters that were used to develop the ANN model and RSM equation.

**Table 2.** Certificates of the first ANN model for the whole dataset.


To investigate the complex influence of *F*15, the new dataset was constructed by selecting data samples with an *F*<sup>15</sup> less than critical value of 28%, which was demonstrated by Tao [32], and its characteristics are listed in Table 3. The new dataset was divided again into three groups with the same portion of each site and with similar statistical factors. Therefore, data from each earthquake contributed to the training, testing, and validating phase. In this step, around 15%, 15%, and 70% of the dataset equated to 27, 27, and 129 samples used for the testing, validating, and training processes, respectively. The characteristics of this new ANN model are presented in Table 3. Also, Table 4 presents the certificate of the second model. It can be seen in Table 4 that the *R* values for all groups of datasets were around 90%.

**Table 3.** Characteristics of database with *F*<sup>15</sup> ≤ 28% used for the second developed ANN model and RSM equation.




#### *5.3. The RSM Equations for Predicting DH*

First, an RSM was conducted to drive an equation on the basis of the first ANN model. Therefore, the DOE introduced by Box et al. [43] for the second-degree polynomial with cross terms was employed to provide 54 coded values to cover the full range of parameters. Through the first developed ANN model, the response value (herein referred to as *DH*) in coded form was calculated. Thereafter, the RSM equation with 28 terms according to the second-degree polynomial with cross terms was derived; however, through hypothesis testing considering the *p*-value, some terms were eliminated. Then, the RSM was applied repeatedly to achieve the final equation. The following equation was consequently developed with 22 terms to correlate the *DH* caused by liquefaction to the six input

parameters for the whole range of parameters in this study, without any limitations on the range of the *F*<sup>15</sup> value:

$$D\_{H} = a\_{0} + a\_{1}M\_{w} + a\_{2}\mathcal{W} + a\_{3}T\_{15} + a\_{4}F\_{15} + a\_{5}(D50\_{15}) + a\_{6}(CAV\_{5}) + a\_{7}M\_{w}^{\;2} + a\_{8}\mathcal{W}^{2} + a\_{9}T\_{15}^{\;2}$$

$$+ a\_{10}(F\_{15})^{2} + a\_{11}(D50\_{15})^{2} + a\_{12}(CAV\_{5})^{2} + a\_{13}M\_{w}T\_{15} + a\_{14}M\_{w}F\_{15} + a\_{15}M\_{w}(D50\_{15}) + a\_{16}\mathcal{W}T\_{15} \tag{13}$$

$$+ a\_{17}\mathcal{W}F\_{15} + a\_{18}\mathcal{W}(D50\_{15}) + a\_{19}\mathcal{T}\_{15}\mathcal{F}\_{15} + a\_{20}\mathcal{T}\_{15}\mathcal{D}50\_{15} + a\_{21}\mathcal{F}\_{15}\mathcal{CAV}\_{5}$$

Characteristics of the RSM equation: R2 = 87.22%, R2 (predicted) = 78.73%, and R2 (adjust) = 83.99%.

The coefficient of determination (R2) illustrates how well the curves fit on the data points. In addition, the R2 (adjust) demonstrates the percentage of variation defined by the independent variables that affect the dependent variable, herein referred to as *DH*. Also, the R<sup>2</sup> (predicted) defines how well a correlation is able to predict the target for a new observation. The values of coefficients *a*<sup>0</sup> to *a*<sup>21</sup> are listed in Table 5.

Coefficient *a***<sup>0</sup>** *a***<sup>1</sup>** *a***<sup>2</sup>** *a***<sup>3</sup>** *a***<sup>4</sup>** *a***<sup>5</sup>** *a***<sup>6</sup>** *a***<sup>7</sup>** Value 0.9174 −1.6737 2.6172 0.7685 −1.0865 −1.8952 1.3425 −0.36369 Coefficient *a***<sup>8</sup>** *a***<sup>9</sup>** *a***<sup>10</sup>** *a***<sup>11</sup>** *a***<sup>12</sup>** *a***<sup>13</sup>** *a***<sup>14</sup>** *a***<sup>15</sup>** Value −0.3733 −0.0678 −0.7474 −0.4060 0.0258 −0.3766 0.2579 −0.59428 Coefficient *a***<sup>16</sup>** *a***<sup>17</sup>** *a***<sup>18</sup>** *a***<sup>19</sup>** *a***<sup>20</sup>** *a***<sup>21</sup>** Value 0.3566 −0.4549 0.4603 −0.6531 0.6011 −0.5063

**Table 5.** Coefficients of Equation (13).

Second, the second RSM equation with the final 21 terms was derived based on the second developed ANN model in this study through the same process as that for the first RSM equation. The second-degree polynomial with cross terms was applied in conjunction with Box and Behnken's DOE on the basis of the second ANN model. Then, a hypothesis test with the same *p*-value was conducted to provide the second RSM equation for *F*<sup>15</sup> less than 28% (Equation (14)). Table 6 presents the coefficients of the second RSM equation.:

$$\begin{aligned} D\_{H} &= a\_{0} + a\_{1}M\_{\overline{w}} + a\_{2}\mathcal{W} + a\_{3}T\_{15} + a\_{4}F\_{15} + a\_{5}(\text{D50}\_{15}) + a\_{6}(\text{CAV}\_{5}) + a\_{7}M\_{\overline{w}}^{2} + a\_{8}\mathcal{W}^{2} + a\_{9}T\_{15}^{2} \\ &+ a\_{10}(F\_{15})^{2} + a\_{11}(\text{D50}\_{15})^{2} + a\_{12}(\text{CAV}\_{5})^{2} + a\_{13}M\_{\overline{w}}\mathcal{W} + a\_{14}M\_{\overline{w}}F\_{15} + a\_{15}M\_{\overline{w}}(\text{CAV}\_{5}) \\ &+ a\_{16}\mathcal{W}(\text{D50}\_{15}) + a\_{17}\mathcal{W}(\text{CAV}\_{5}) + a\_{18}\begin{aligned} F\_{15} &\ F\_{15} + a\_{19}\mathcal{T}\_{15} \left(D50\_{15}\right) + a\_{20}\mathcal{T}\_{15} \left(D50\_{15}\right) \end{aligned} \end{aligned} \tag{14}$$

Characteristics of the second RSM equation: R<sup>2</sup> = 88.51%, R2 (predicted) = 50.95%, and R<sup>2</sup> (adjust) = 78.09%.


**Table 6.** Coefficients of Equation (14).

It must be noted that to use the derived RSM equation, the main values of input parameters must be exchanged with coded values by the function presented in Equation (15). Those coded values must then be put into the RSM equation to achieve the results. In other words, the RSM equation declares the relationship between the coded value of parameters and responses (output); coded values have a range from −1 to 1.

$$\text{Coded value} = \frac{\text{Real value} - \text{mean value}}{\text{Max value} - \text{Min value}} \times \text{2} \tag{15}$$

The max, min, and mean values of the parameters are listed in Tables 1 and 3 for both RSM equations presented in this study.

#### **6. Comparison of RSM Equations with Extra Models**

Chu et al. [47] analyzed five liquefied sites during the Chi-Chi-1999 earthquake in Taiwan, all in the near-fault region from five sites in two districts of Wufeng and Natu as they are illustrated in Figures 1 and 2. Table 7 presents these parameters' values from the sites. Based on these data samples, the predicted results of the RSM equations are compared to Youd et al. [20], Javadi et al. [27], and Rezania et al. [36]. In total, 28 sites (for which the necessary parameters for the ANN model and the RSM equation are reported by Chu et al. [45]) are illustrated in Table 8. The sample numbers from 1 to 26 are from Wufeng's sites and samples of 27 and 28 are belong to Nantu's site. The capability and accuracy of the first RSM equation was demonstrated using all 28 site samples from the Chi-Chi earthquake, including samples with *F*<sup>15</sup> from 13% to 48.5%.

**Table 7.** Model parameters and measured *DH* at sites affected by the 1999 Chi-Chi earthquake.


The results of comparison between the predicted values and measured values are summarized in terms of the root mean square error (*RMSE*), mean absolute error (*MAE*), and *R* in Table 8. It is clear that the larger *R* and smaller *RMSE* and *MAE* reveal higher accuracy of predicted results.

$$RMSE = \sqrt{\frac{\sum\_{N} \left(X\_m - X\_P\right)^2}{N}} \tag{16}$$

$$MAE = \frac{\sum\_{N} |X\_m - X\_P|}{N} \tag{17}$$

where *N* is the number of samples, *Xm* is the measured value, and *XP* is the predicted value.


**Table 8.** Performance certificate of first RSM equation in comparison with extra available models.

Furthermore, all samples with an *F*<sup>15</sup> greater than 28% were eliminated from the Chi-Chi earthquake cases, and 16 samples consequently remained (samples number 1 to 14 as well as numbers 27 and 28, as can be seen in Table 7). Then, the second RSM equation was validated by applying it to these samples in comparison with the extra three models. Table 9 summarizes the results of all models.

**Table 9.** Performance certificate of second RSM equation, on the basis of samples with *F*<sup>15</sup> ≤ 28% in comparison with extra available models.


Figures 6 and 7 visualize the comparison between both RSM equations developed in the present study and three extra models with measured data from sites of the Chi-Chi earthquake. Twenty-eight data points are evaluated in Figure 6 for whole range of the parameters. Meanwhile, Figure 7 illustrates the comparison for data points with *F*<sup>15</sup> values of less than 28% at 16 data points.

**Figure 6.** Comparison between the first RSM equation and three extra models with 28 data points measured from sites of Chi-Chi earthquake.

**Figure 7.** Comparison between both RSM equations and three extra models with 16 data points including *F*<sup>15</sup> of less than 28% measured from sites of the Chi-Chi earthquake.

#### **7. Results and Discussion**

The previous sections have compared the first RSM model, which belongs to the full range of parameters, and the second RSM model, which was derived for samples whose *F*<sup>15</sup> values were less than 28%, with three extra well-known models. The models were examined using new data from the Chi-Chi earthquake, which were not included in the two datasets to establish the two RSM models. As can be seen in Table 8, the RSM equation of the whole range of parameters indicated a higher *R* value of 0.683, in comparison with the extra models whose values were 0.433, −0.74, and 0.514. Furthermore, the RSM model comprises lower *MAE* and *RSME* values of 0.3 and 0.37, respectively, compared to 0.49 and 0.7 for Rezania et al., 1.04 and 1.19 for Javadi et al., and 3.77 and 4.37 for Youd et al. Therefore, among all of them, the RSM model provided prediction with higher accuracy.

On the other hand, as can be seen in Table 9, by considering samples with *F*<sup>15</sup> less than 28%, the model of Youd et al. provided the highest *R* value of 0.934, closely followed by the second and the first RSM models with *R* values equal to 0.891 and 0.846 respectively. Table 9 also illustrates the MAE and RSME criteria values for all models for samples with a limited value for *F*<sup>15</sup> less than 28%. The values of the *MAE* and *RSME* in the second RSM model—0.29 and 0.39, respectively—indicates the highest accuracy and performance in comparison with the others. In addition, the model of Rezania et al., with 0.42 and 0.57, illustrated lower values for the *MAE* and *RSME*, respectively. Further, Javadi et al. with 1.2 and 1.3, and Youd et al. with 4.84 and 5.34 provide less accuracy for predicting *DH*.

The comparison between the two models developed in this study and the extra three models demonstrates that the second RSM model provided a reasonable correlation and the lowest error. The results indicated that the RSM is a highly efficient tool to perform a liquefaction hazard analysis. Furthermore, performance of the model is increased by taking into account the complex influence of *Fc* by eliminating an *F*<sup>15</sup> larger than 28% and even by decreasing the number of samples in the dataset.

Another major advantage of the presented models is their consideration of earthquake aspects, such as the near-fault zone, the frequency of earthquake motion, and the causative fault type, by estimating and adding the *CAV*<sup>5</sup> parameter to the dataset. As can be seen from Figures 6 and 7, among all models that were considered in the present study to calculate *DH* without any limitation on the parameters' value, the model of Youd et al. was overpredicted. Meanwhile, Youd et al.'s model provided poor and overpredicted results for samples with a limited value of *F*<sup>15</sup> less than 28%. Additionally, considering samples with a limited *F*<sup>15</sup> value shows the first RSM and the model of Javadi et al. present an overpredicted value for *DH*. Furthermore, second RSM equation and model of Rezania et al. underpredicted *DH* in their predictions.

There are some limitations for applying both first and second RSM equations as follow:


#### **8. Summary and Conclusions**

The determination of lateral displacement due to liquefaction caused by an earthquake (*DH*) is the most important aspect of liquefaction hazard analysis. There are two main types of conditions according to the topography of the sites: free-face and sloping ground conditions. First, the parameter of corrected absolute velocity (*CAV*5) of sites was calculated due to it being the most efficient and sufficient parameter for the assessment of liquefaction caused by earthquakes [21], and it was added to develop the dataset to cover all aspects of earthquakes, including the frequency content of earthquake motions and the causative fault type of earthquakes. Then, a statistical parametric analysis was performed by estimating the correlation coefficient (R) between all input parameters and output as *DH*. To achieve a more capable and accurate model, based on the estimated values for *R*, the horizontal distance from a site to the seismic energy source (*r*) and ground slope (*S*) was eliminated from the original dataset due to poor correlations to the target. Therefore, the final dataset was created for free-face condition sites.

The significant aspects of earthquakes, such as the near-fault region, frequency content, and causative fault type of earthquakes, which are included in the model established by Kramer et al. [21], were considered by taking *CAV*<sup>5</sup> into account. To investigate the complex effect of fine content, the main dataset was divided into two subsets. The first dataset included the whole range of parameters, and in the second one, all samples with average fine content in the liquefiable layer (*F*15) larger than 28% were removed from the dataset, in line with Tao [32]. Furthermore, the RSM was applied to develop two equations in order to the first and the second dataset to examine its performance to assess liquefaction. In the end, the two presented models in this study were compared to three available models to

demonstrate their capability and accuracy with regard to predicting *DH* in free-face conditions in a near fault zone case history of the Chi-Chi earthquake.

The present study highlights the importance of earthquake aspects, especially *CAV*<sup>5</sup> as the most sufficient and efficient intensity to liquefaction hazard assessments. In addition, the RSM is a strong tool for the evaluation of complex non-linear phenomena such as liquefaction.

The results also confirm the complicated influence of *F*<sup>15</sup> on the whole range, and they provide significant enhancements to the performance of the model by considering samples with an *F*<sup>15</sup> less than 28% as a critical value defined by Tao [32]. One of the most remarkable results, which shows the complex influence of fine content on evaluation of *DH*, is that the second model demonstrated higher accuracy and capability, even though it was developed using a database with fewer samples than the first model.

**Author Contributions:** Conceptualization, X.T.; methodology; Q.Y. and X.T.; software, N.P.; validation, N.P.; formal analysis, N.P.; investigation, Q.Y. and X.T.; resources, N.P.; data curation, N.P.; writing—original draft preparation, N.P.; writing—review and editing, N.P.; visualization, N.P.; supervision, Q.Y. and X.T.; project administration, X.T.; funding acquisition, Q.Y.

**Funding:** This research was funded by National Natural Sciences Foundation of China, grant number 51639002 and National Key Research & Development Plan, grant number 2018YFC1505305 and "The APC was funded by State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China.

**Acknowledgments:** The authors are grateful for the technical and financial support Provided by National Natural Sciences Foundation of China Granted No. 51639002 and National Key Research & Development Plan under Grant No. 2018YFC1505305. State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Simulations of Wave-Induced Soil Erosion in Silty Sand Seabeds**

#### **Zhen Guo, Wenjie Zhou, Congbo Zhu, Feng Yuan \* and Shengjie Rui**

College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China; nehzoug@163.com (Z.G.); zhouwenjiesd@163.com (W.Z.); 21812141@zju.edu.cn (C.Z.); ruishengjie@zju.edu.cn (S.R.)

**\*** Correspondence: yuanfen5742@163.com; Tel.: +86-137-3547-4967

Received: 10 January 2019; Accepted: 18 February 2019; Published: 20 February 2019

**Abstract:** Silty sand is a kind of typical marine sediment that is widely distributed in the offshore areas of East China. It has been found that under continuous actions of wave pressure, a mass of fine particles will gradually rise up to the surface of silty sand seabeds, i.e., the phenomenon called wave-induced soil erosion. This is thought to be due to the seepage flow caused by the pore-pressure accumulation within the seabed. In this paper, a kind of three-phase soil model (soil skeleton, pore fluid, and fluidized soil particles) is established to simulate the process of wave-induced soil erosion. In the simulations, the analytical solution for wave-induced pore-pressure accumulation was used, and Darcy flow law, mass conservation, and generation equations were coupled. Then, the time characteristics of wave-induced soil erosion in the seabed were studied, especially for the effects of wave height, wave period, and critical concentration of fluidized particles. It can be concluded that the most significant soil erosion under wave actions appears at the shallow seabed. With the increases of wave height and critical concentration of fluidized particles, the soil erosion rate and erosion degree increase obviously, and there exists a particular wave period that will lead to the most severe and the fastest rate of soil erosion in the seabed.

**Keywords:** wave action; silty sand; seepage flow; soil erosion; pore-pressure accumulation; three-phase soil model

#### **1. Introduction**

Silty sand is widely distributed in the eastern coast of China, among which the most representative area is the Yellow River subaqueous delta. According to the in-situ survey data [1], the silty sand sediment (typical median particle size less than 50.00 μm, silt content over 80%) accounts for 90% of the northeast of the delta. There commonly exists a kind of hard crust with a thickness of 2.00–3.00 m in the shallow stratum of seabeds. Sumer et al. [2] presented the results of an experimental investigation of the complete sequence of sediment behavior beneath progressive waves and reported a similar hard crust in sandy seabeds. The main reason for the formation of hard crust was thought to be the compaction or solidification of sand layers induced by waves. However, for silty sand seabeds, the coarse and fine particles coexist and the particle size distribution varies greatly. Thus, the inner mechanism becomes different and complicated. Under wave actions, fine particles filling in the pore space tend to move with the seepage flow, but the coarse particles remain in their initial positions. This characteristic has been verified by the previous work of Shi [3]. Using a scanning electron microscope, Shi [3] investigated the micro-structures of the hard crust, and found that the hard crust is constituted of uniform coarse particles and a few fine particles. Based on a lot of field and experimental tests, Jia et al. [1] pointed out that the hard crust is mainly caused by the wave-induced reformation and erosion of the sediments near the surface. As shown in Figure 1a [4], under continuous wave actions, plumes of sediment deposit over the seabed surface due to the upward movement of fine particles. Figure 1b [4] shows the micro-holes in the silty sediment as the result of fine particle transportation. This phenomenon is also named "seabed coarsening". The seabed coarsening phenomenon commonly appears in shallow seabeds, but currently suitable theoretical or numerical models are still lacking for the wave-induced erosion process of silty sand seabeds. The coarsening phenomenon of the seabed will lead to the increase of soil permeability, which is the most important effect that can significantly affect the potential and the depth of seabed liquefaction. In addition, the mechanical properties of seabed soil will also be changed when seabed coarsening is occurred.

(**a**) (**b**)

**Figure 1.** Plumes of sediment and micro-holes in silty sediment seabed: (**a**) The plumes of sediment on the silty sediment surface; (**b**) The micro-holes due to erosion. [4].

As mentioned above, the soil erosion is induced by the seepage flow within the seabed under wave actions. Under the extreme wave condition, the excess pore-pressure is always large enough for the occurrence of soil liquefaction, and thus soil particles will be repositioned and reconsolidated [5]. For the normal wave condition, the wave height is small and continuous seepage flow can be induced, so some fine particles in the silty sand seabed move upwards under the seepage force, and the coarsening phenomenon will emerged in shallow seabeds [6], as shown in Figure 2. It is also pointed out that the hydrodynamic condition plays a significant role in topography construction and seabed erosion process.

**Figure 2.** Two mechanisms of the wave-induced pore-pressure and the erosion process.

In this paper, a three-phase soil model (soil skeleton, pore fluid, and fluidized soil particles) was established to study the soil erosion process induced by waves in the silty sand seabed. In the numerical simulation, the Darcy flow law, mass conservation, and generation equations were coupled into COMSOL Multiphysics [7] to perform the studies. COMSOL Multiphysics is a kind of finite element method (FEM) software which is developed by COMSOL INC found in Stockholm, Sweden. Jeng et al. [8] discussed two mechanisms for wave-induced pore pressures in a porous seabed, i.e., oscillatory, residual excess pore pressures, and an analytical solution for the wave-induced residual pore pressure was derived. Using the residual pore-pressure analytical solution [8], the process of wave-induced soil erosion was investigated. Then, the parametric studies were performed to study the influences of wave height, wave period, and critical concentration of fluidized particles on the erosion process in the seabed. It is found that the most significant soil erosion mainly occurred at the shallow seabed. With the increases of wave height and critical concentration of fluidized particles, the soil erosion rate and erosion degree increase obviously, and there exists a particular wave period that will lead to the most severe and the fastest rate of soil erosion in the seabed.

#### **2. Analytical Solution for Wave-Induced Pore-Pressure Accumulation**

Generally speaking, based on the generation mechanism, as shown in Figure 2, the total excess pore-pressure is composed of the oscillatory pore-pressure and the residual pore-pressure when waves propagate along the seabed surface [5,9–12], and it can be expressed by

$$P = P\_{\text{osc}} + P\_{\text{res}} \tag{1}$$

where *P*osc is the oscillatory pore-pressure corresponding to the elastic deformation of the soil skeleton. *P*osc fluctuates in both temporal and spatial domains, and the fluctuation is accompanied by the attenuation of the amplitude and phase lag under wave actions [13–15]. *P*res is the residual pore-pressure that is period-averaged, and is the result of accumulated plastic deformation of the soil skeleton. It has been acknowledged recently that with the accumulation of pore-pressure, continuous seepage flow appears near the seabed surface and may lead to obvious particle migration [16–18].

Many studies have been performed for the accumulation of excess pore-pressure in the seabed induced by waves [8,19–22]. According to Jeng et al. [8], for the waves, according to linear wave theory, the residual pore-pressure in infinite thickness seabed can be derived based on Biot's consolidation equation in one-dimension [23], and the analytical solution can be expressed as

$$P\_{\rm res} = \frac{2A}{c\_v \lambda\_k^{-3}} [1 - (\frac{\lambda\_k z}{2} + 1) \exp(-\lambda\_k z) - \frac{1}{\pi} \int\_0^\infty \frac{\exp(-r\varepsilon\_v \lambda\_k^{-2} t)}{r(r+1)^2} \sin(\sqrt{r} \lambda\_k z) dr] \tag{2}$$

$$A = \frac{\gamma t (1 + 2K\_0)}{3T} \left[ \frac{3P\_b k\_s}{\beta (1 + 2K\_0) \gamma \prime} \right]^{1/\eta} \tag{3}$$

$$
\lambda\_k = \frac{k\_s}{\eta} \tag{4}
$$

where *cv* is the consolidation coefficient, *K0* is the coefficient of lateral earth pressure, *β* and *η* are empirical constants, which can be confirmed based on the soil type and the relative density [24], *ks* is the wave number, *T* is the wave period, and *Pb* is the amplitude of the dynamic wave pressure on the seabed surface.

#### **3. Theoretical Model for Soil Erosion Process**

#### *3.1. Definition of Three-Phase Soil Model*

As shown in Figure 3, under the effect of wave-induced seepage, the transportation of the fine particles will be induced. In this paper, a kind of three-phase soil model is defined and used to simulate the transportation process of fine particles. The three-phase model was first proposed by Vardoulakis et al. [25] to analyze the sand production problem. Accordingly, the soil element is defined to be the combination of the soil skeleton (*s*), pore fluid (*f*), and fluidized soil particles (*fs*), which can be expressed as

$$d\mathcal{W} = d\mathcal{W}\_f + d\mathcal{W}\_{fs} + d\mathcal{W}\_s \tag{5}$$

where *dW*, *dWf*, *dWf s*, and *dWs* are the volumes of the soil element, soil skeleton, pore fluid, and fluidized soil particles, respectively. The masses of soil element, soil skeleton, pore fluid, and fluidized soil particles are represented by *dM*, *dMf*, *dMfs*, and *dMs*, respectively.

**Figure 3.** Three-phase theoretical model for the seabed soil.

In the three-phase model, the velocities of the three phases are

$$
v\_{fs} = v\_f = \overline{v} \tag{6}$$

$$v\_{\mathbb{S}} = 0 \tag{7}$$

where *vf s*, *vf* , *v*, *vs* are the velocities of the fluidized soil particles, pore fluid, the mixture (pore fluid and fluidized soil particles), and soil skeleton, respectively.

The concentration of the fluidized soil particles *c* can be expressed by

$$\mathcal{L} = \frac{d\mathcal{W}\_{fs}}{d\mathcal{W}\_{fs} + d\mathcal{W}\_f} \tag{8}$$

The soil porosity *ϕ* can be defined as

$$\varphi = \frac{d\mathcal{W}\_f + d\mathcal{W}\_{fs}}{d\mathcal{W}} \tag{9}$$

The density of the mixture *ρ* is

$$
\overline{\rho} = (1 - c)\rho\_f + c\rho\_s \tag{10}
$$

where *ρf, ρ<sup>s</sup>* are the densities of the pore fluid and the solid skeleton.

The apparent density of the fluidized soil particles can be defined as

$$
\overline{\rho}\_{fs} = \frac{dM\_{fs}}{d\mathcal{W}} = c\rho\rho\_s \tag{11}
$$

The volume discharge rate *q* and the velocity *v* of the mixture are

$$
\overline{q} = \frac{d\overline{W}}{dSdt} \tag{12}
$$

$$\overline{v} = \frac{d\overline{W}}{d\overline{S}dt} = \frac{d\overline{W}}{\overline{q}dSdt} = \frac{\overline{q}}{\overline{q}}\tag{13}$$

where *dW* is the volume of the mixture through the cross-sectional *dS* within *dt* time, *dS* is the pore part of *dS*.

*J. Mar. Sci. Eng.* **2019**, *7*, 52

#### *3.2. Mass Conservation Equations*

Vardoulakis et al. [25] and Sterpi [26] introduced the mass conservation equation of the three-phase in one-dimension shown as

$$\frac{\partial \rho\_{\mathfrak{a}}}{\partial t} + \frac{\partial}{\partial z} (\rho\_{\mathfrak{a}} v\_{\mathfrak{a}}) = \dot{m}\_{\mathfrak{a}} \tag{14}$$

where . *m<sup>α</sup>* is the mass generation term, which means the mass generation rate of phase *α* (the phase α can represent the fluidized particles phase with subscript *fs* or solid phase with subscript *s* or fluid phase with subscript *f*), and *∂ρα <sup>∂</sup><sup>t</sup>* is the density change rate with time of phase *α*.

In detail, the three phases can be expressed as follows and the related diagrams are shown in Figure 4.

**Figure 4.** Mass conservation of the three phases: (**a**) Fluidized particles; (**b**) soil skeleton; (**c**) Pore fluid.

#### (1) Fluidized soil particles

Combining Equation (6), (11), and (14), the mass conservation equation of the fluidized soil particles can be expressed as

$$\frac{\partial(c\,\rho)}{\partial t} + \frac{\partial}{\partial z}(c\,\rho\overline{v}) = \frac{\dot{m}}{\rho\_s} \tag{15}$$

where . *<sup>m</sup>* <sup>=</sup> . *mer* <sup>−</sup> . *mdep*, . *mer* is the rate of eroded mass and . *mdep* is the rate of deposited mass.

#### (2) Soil skeleton

Here, we divided the soil element into the solids (index 1) and the mixture (index 2). According to Equation (14), the mass conservation equations of the two phases are

$$-\dot{m} = \frac{\partial \rho\_1}{\partial t} + \frac{\partial}{\partial z}(\rho\_1 v\_1) \tag{16}$$

$$
\dot{m} = \frac{\partial \rho\_2}{\partial t} + \frac{\partial}{\partial z}(\rho\_2 v\_2) \tag{17}
$$

where *ρ1, ρ<sup>2</sup>* are the densities of the soil phase and the mixture phase, *v1, v2* are the velocities of the soil phase and the mixture phase.

The mass conservation equation of the soil skeleton is

$$\frac{\partial \varrho}{\partial t} = \frac{\dot{m}}{\rho\_s} = \frac{\dot{m}\_{er} - \dot{m}\_{dep}}{\rho\_s} \tag{18}$$

Using Equation (18), Equation (15) can be re-expressed by

$$\frac{\partial \overline{\phi}}{\partial t} = \frac{\partial (c\overline{\phi})}{\partial t} + \frac{\partial}{\partial z}(c\overline{\phi}\overline{v}) \tag{19}$$

#### (3) Pore fluid

Combining Equation (10) and (15), Equation (17) can be transformed into

$$
\frac{
\partial
}{
\partial t
}[(1-c)\varphi] + \frac{
\partial
}{
\partial z
}[\varphi(1-c)\overline{\varphi}] = 0
\tag{20}
$$

With Equation (19), Equation (20) can be re-expressed by

$$\frac{\partial(\varrho\overline{\upsilon})}{\partial z} = 0 \tag{21}$$

Thus, the simplifications of these three equations are

$$\begin{cases} \begin{array}{c} \frac{\partial \varphi}{\partial t} = \frac{\dot{m}}{\rho\_\*} \\\\ \frac{\partial \varphi}{\partial t} = \frac{\partial (c\varphi)}{\partial t} + \frac{\partial}{\partial \overline{z}} (c\varphi \overline{v}) \\\\ \frac{\partial (\varphi \overline{v})}{\partial \overline{z}} = 0 \end{array} \tag{22}$$

There are four basic variables (*ϕ*, . *<sup>m</sup>*, *<sup>c</sup>*, *<sup>v</sup>*) in Equation (22), and a constituted equation for . *m* is needed to solve the problem.

#### *3.3. Constitutive Laws of Mass Generation*

The rate of the soil erosion . *mer* can be expressed by

$$
\dot{m}\_{cr} = \rho\_s \lambda (1 - \varrho) \mathbf{c} ||\overline{\eta}|| \tag{23}
$$

where *λ* is the parameter used to describe the spatial frequency of the potential erosion starter points in the soil skeleton of the porous medium and can be obtained using experiments [25]. It can be seen that . *mer* is proportional to *c*, which means the erosion process can go on until *c* is equal to 0. The particle deposition takes place in parallel with the particle erosion. According to Vardoulakis et al. [25], the particle deposition rate can be expressed by

$$
\dot{m}\_{dep} = \rho\_s \lambda (1 - \varphi) \frac{c^2}{c\_{\mathcal{T}}} ||\overline{\eta}|| \tag{24}
$$

Combining Equations (23) and (24), the net particle erosion . *m* can be expressed by

$$
\dot{m} = \dot{m}\_{cr} - \dot{m}\_{dep} = \rho\_s \lambda (1 - \varrho) \left( \varepsilon - \frac{c^2}{c\_{cr}} \right) \|\overline{q}\|\tag{25}
$$

#### *3.4. Darcy Flow Law*

With the loss of fine particles in the erosion process, the grain size distribution of the silty sand will be changed and the soil porosity will be increased. Grain size distribution of sand affects its permeability. It is known that poorly-graded soil has higher porosity and its permeability is larger than that of the well-graded soil, in which smaller grains tend to fill the voids between larger grains. According to the Carman-Kozeny equation [27], the relationship between the soil permeability and the porosity can be described as

$$k = K \frac{\varrho^3}{\left(1 - \varrho\right)^2} \tag{26}$$

where *k* is the soil permeability, *K* is the reference permeability.

The seepage flow under hydraulic gradient can be described by Darcy flow law [28], shown as

$$\overline{\eta} = -\frac{k}{\eta\_k \overline{\rho}} \cdot \frac{\partial P}{\partial z} \tag{27}$$

where *η<sup>k</sup>* is the kinematic viscosity of the mixture of pore fluid and fluidized particles.

#### *3.5. Governing Equations for Soil Erosion*

By including mass conservation equations, mass generation law, and Darcy flow law, the governing equations for the soil erosion process induced by waves in one-dimension are shown in Equation (28).

$$\begin{cases} \begin{aligned} \frac{\partial \overline{\rho}}{\partial t} &= \frac{\partial (c\overline{\rho})}{\partial t} + \frac{\partial (c\overline{\rho})}{\partial z} \\\\ \frac{\partial \overline{\rho}}{\partial t} &= \lambda (1 - \overline{\rho}) (c - \frac{c^2}{\zeta\_{cr}}) ||\overline{\eta}|| \end{aligned} \end{cases} \tag{28}$$
 
$$\begin{aligned} \frac{\partial (\overline{\eta})}{\partial z} &= 0 \\\\ \overline{\eta} &= -\frac{K\rho^3}{(1 - \rho)^2 \eta\_k \overline{\rho}} \cdot \frac{\partial P}{\partial z} \end{aligned} \tag{29}$$

In Equation (28), the basic variables are only *ϕ*, *c*, *P*, and all of which are the functions of time *t* and position *z*.

#### **4. Numerical Implement of Seabed Erosion Model and Simulations**

In this section, a numerical model was established to analyze the erosion process of silty sand seabeds induced by waves. COMSOL Multiphysics is a kind of general-purpose simulation software for FEM modelling in all fields of engineering and scientific research [7]. In this paper, the Partial Differential Equation (PDE) module was used for the secondary development. In detail, the numerical implement process can be described as follows. Firstly, the residual pore-pressure in the seabed induced by waves can be obtained using Equation (2) proposed by Jeng [8]. The distribution of the residual pore-pressure is inputted into the seabed erosion model. Then, with full drainage conditions on the seabed surface and the impermeable seabed bottom, the Darcy seepage process can be solved. For the seabed erosion model, the PDE module in COMSOL is used to solve Equation (28), thus the erosion process (changes of *ϕ, c*) can be obtained. In the numerical model, the Lagrange shape function and the quadratic element order were adopted. The backward difference method was selected to discretize the time domain and the Newton-Raphson method was used to solve the governing equations iteratively. To satisfy the request of convergence, the time step Δ*t* satisfy

$$
\Delta t \le \frac{l}{\sqrt{E/\rho\_s}}\tag{29}
$$

where *l* is length of the minimum element, *E* is elastic modulus of soil.

In the numerical model, the geometry of seabed depth *ds* is equal to 30.00 m and the average mesh size is 0.1 m. More parameters can be listed as follows: water depth *dw* = 10.00m, wave height *H* = 2.00 m, wave period *T* = 5.00 s, wave length *L* = 36.59 m. According to the judgement criterion about the seabed depth [12], *ds*/*L* = 0.82>0.3, and thus the depth of seabed can be treated as infinite thickness. For the soil condition, the initial porosity *ϕ<sup>0</sup>* = 0.42, initial concentration of the fluid soil particles *c0* = 0.001. More details can be found in Table 1. For a typical wave condition, a series of numerical studies have been performed. It is known that the wave-induced erosion is not only associated with soil properties, but also closely related to wave characteristics. So, the influences of wave height *H*, wave period *T*, and critical concentration of the fluidized soil particles *ccr* on the process of wave-induced erosion were discussed. The simulation cases are listed in Table 2.


**Table 1.** Parameters used in the numerical model for the typical wave condition case.

**Table 2.** Calculation cases of the parametric analyses.


#### **5. Time Characteristics of Wave-Induced Soil Erosion Process**

To investigate the time characteristics of the wave-induced soil erosion process, a typical wave condition case under normal sea state was analyzed in the simulation. The wave acting time *t* was selected as 1 h, 2 h, 5 h,10 h, 24 h, 2 d, 5 d, 15 d, 30 d (h is one hour and d refers to one day), respectively. The distributions of the oscillatory pore-pressure and the residual pore-pressure in the seabed are shown in Figure 5. As shown in Figure 5a, the dimensionless maximum oscillatory pore-pressure |*Posc*|/*Pb* decreases from 1.00 on the seabed surface to 0 at the −30.00 m depth. The liquefaction of the seabed can be divided into the oscillatory and residual liquefactions [5]. According to Jeng et al. [8] and Okusa [29], the criterions of oscillatory and residual liquefactions are *Posc σ* 0 <sup>≥</sup> 1 and *Pres σ* 0 ≥ 1, respectively (*σ* <sup>0</sup> is the effective vertical stress of soil). Figure 5b indicates that the oscillatory liquefaction will not occur under the typical wave condition. Figure 5c shows the evolution of the residual pore-pressure along depths. It is noted that the residual pore-pressure develops gradually with the extension of wave acting time and tends to be stable. The maximum value of *Pres* occurs at about −5 m to −10 m (below the seabed surface) depth in the whole process of wave actions. In Figure 5d, it also reveals that there is no potential soil liquefaction in the seabed with the accumulation of *Pres*. Under normal sea state, the soil erosion is the common behavior for the silty sand seabed.

In the erosion process, part of the soil skeleton is transformed into fluidized particles, which remain in suspension under the effect of seepage flow, and thus the concentration of fluidized particles will be increased. Figure 6 shows the variations of *c* along depth for different wave acting times. It shows that the maximum value of *c* occurs on the seabed surface in the erosion process. For the shallow depth (within −2.00 m), *c* increases from the initial value 0.001 to the critical value 0.30 and then keeps a stable state. When the wave acting time *t* = 30 d, the seabed depth affected by wave-induced erosion is up to −4.00 m.

**Figure 5.** Distributions of the oscillatory pore-pressure and the residual pore-pressure: (**a**) Vertical distribution of |*Posc*|/*Pb*; (**b**) vertical distribution of (*Pb* − |*Posc*|)/*σ* <sup>0</sup>; (**c**) vertical distribution of *Pres* for different times; (**d**) vertical distribution of *Pres*/*σ* <sup>0</sup> for different times.

**Figure 6.** Variations of the concentration of moving particles with the increase of wave acting time: (**a**) Diagram of three-dimensions; (**b**) diagram of two-dimensions.

The soil porosity increases with the loss of fine particles during the erosion process. It can be seen in Figure 7 that the soil porosity gradually increases with the extension of wave acting time at shallow depths (within −5.00 m), and the soil porosity in deep depths remains almost constant. When the wave acting time is less than 24 d, the maximum value of soil porosity occurs on the seabed surface. After 24 d, the most severe erosion occurs at the depth of about −0.50 m, and the soil porosity keeps the value of 0.55 on the seabed surface. It illustrates that the greatest loss of fine particles occurs at approximately −0.50 m depth. The evolution of soil permeability in the erosion process is shown in Figure 8. It is shown that *k*/*k0* increases with the extension of wave acting time and the maximum value reaches 4.10 at −0.50 m depth after 24 d. When *t* = 30 d, the depth with *k*/*k0* over 2.00 is around −2.30 m. These results indicate that the soil permeability increases significantly with the extension of wave acting time at the shallow seabed.

**Figure 7.** Variations of the soil porosity with the increase of wave acting time: (**a**) Diagram of three-dimensions; (**b**) diagram of two-dimensions.

**Figure 8.** Variations of the soil permeability with the increase of wave acting time: (**a**) Diagram of three-dimensions; (**b**) diagram of two-dimensions.

Two physical quantities *<sup>∂</sup><sup>c</sup> <sup>∂</sup><sup>t</sup>* and *∂ϕ <sup>∂</sup><sup>t</sup>* are introduced in this paper to describe the rate of the soil erosion at every moment, as shown in Figure 9. It can be seen that the erosion rate firstly increases until reaching the peak value, and then gradually decreases. The deeper the seabed soil, the later the peak values of *<sup>∂</sup><sup>c</sup> <sup>∂</sup><sup>t</sup>* , *∂ϕ <sup>∂</sup><sup>t</sup>* can reach and the smaller the peak values of *<sup>∂</sup><sup>c</sup> <sup>∂</sup><sup>t</sup>* , *∂ϕ <sup>∂</sup><sup>t</sup>* . On the seabed surface, the erosion rate reaches the peak value fastest and decreases to negative values, which indicates that the deposition effects play an obvious role in the later stage of the erosion process.

**Figure 9.** Variations of *<sup>∂</sup><sup>c</sup> <sup>∂</sup><sup>t</sup>* , *∂ϕ <sup>∂</sup><sup>t</sup>* at different depths with the increase of wave acting time: (**a**) *<sup>∂</sup><sup>c</sup> <sup>∂</sup><sup>t</sup>* − *t*; (**b**) *∂ϕ <sup>∂</sup><sup>t</sup>* − *t*.

#### **6. Results for Affecting Factors and Interpretations**

#### *6.1. Effect of Wave Height*

Wave height is one of the most important wave parameters, as it directly affects the wave pressure and energy inputted into the seabed [12]. To assess the effect of wave height on the erosion process, the wave height was selected as 1.50 m, 1.75 m, 2.00 m, 2.25 m, and 2.50 m, respectively. Figure 10a shows the distributions of |*Posc*|/*Pb* along depth for different wave heights, which reveal |*Posc*|/*Pb* is only related to soil depth and has no relationship with wave height. In Figure 10b–d, it can be seen that (*Pb* − |*Posc*|)/*σ* <sup>0</sup>, *Pres*, and *Pres*/*σ*<sup>0</sup> increase obviously with the growth of wave height. Compared with the oscillatory pore-pressure, the residual pore-pressure increases more rapidly with the increase of wave height when *t* = 30 d. No oscillatory liquefaction occurs, and the residual liquefaction only occurs with *H* = 2.50 m and *t* = 30 d at shallow depths (within −1.80 m).

Figure 11 shows the evolution of soil porosity with wave height when *t* = 30 d. It can be seen that the soil porosity increases significantly at shallow seabeds with the growth of wave height. The affected depth increases from −2.00 m to −6.00 m when the wave height increases from 1.50 m to 2.50 m. It is also noted that when the wave height is bigger than 2.00 m, the soil erosion on the seabed surface develops rapidly. The effect of wave height on the erosion rate *∂ϕ <sup>∂</sup><sup>t</sup>* is shown in Figure 12. The soil erosion rate at shallow depths increases obviously with the growth of wave height. Similar to Figure 9b, when *H* equals 2.00 m, 2.25 m, and 2.50 m, a negative value of *∂ϕ <sup>∂</sup><sup>t</sup>* appears on the seabed surface at a certain time and then the value becomes positive later in the erosion process. It illustrates that the erosion effect plays a main role again after the deposition effect takes the lead.

**Figure 10.** Distributions of the oscillatory pore-pressure and the residual pore-pressure for different wave heights: (**a**) Vertical distribution of |*Posc*|/*Pb* for different *H*; (**b**) vertical distribution of (*Pb* − |*Posc*|)/*σ* <sup>0</sup> for different *H*; (**c**) vertical distribution of *Pres* for different *H*; (**d**) vertical distribution of *Pres*/*σ* <sup>0</sup> for different *H*.

**Figure 11.** Variations of the soil porosity for different wave heights.

**Figure 12.** Variations of *∂ϕ <sup>∂</sup><sup>t</sup>* for different wave heights: (**a**) *H* = 1.50 m; (**b**) *H* = 2.00 m; (**c**) *H* = 2.25 m; (**d**) *H* = 2.50 m.

#### *6.2. Effect of Wave Period*

The wave length is always related to the wave period and water depth [12]. In this section, the effect of wave period on the wave-induced erosion was studied. The wave period was selected as 2 s, 5 s, 10 s, 15 s, and 20 s. The responses of the pore-pressure for different wave periods are plotted in Figure 13. Figure 13a clearly shows that the maximum value of oscillatory pore-pressure and the affected depth increases significantly with the extension of the wave period. The responses of (*Pb* − |*Posc*|)/*σ* <sup>0</sup>, *Pres* and *Pres*/*σ*<sup>0</sup> for different wave periods are shown in Figure 13b–d, respectively. These three physical quantities first increase and then decrease. There is a competition mechanism between the accumulation and the dissipation of the residual pore-pressure. For the waves with bigger periods, the dissipation of residual pore-pressure becomes relatively obvious. Therefore, there exists a particular wave period corresponding to the maximum residual pore-pressure. The oscillation liquefaction will not occur due to (*Pb* − |*Posc*|)/*σ* <sup>0</sup> always being less than 0.1, and the residual liquefaction appears on the seabed surface when wave period *T*=10s.

The wave period has obvious effects on the soil porosity (Figure 14). The affected depth increases greatly when *T*>5 s. For the soil within the affected depth in the seabed, its porosity increases with the extension of wave period first, and then shows a decreasing trend. The soil erosion is not obvious with a small or big wave period and there exists a particular wave period to make the soil erosion most severe in the seabed. For wave period *T* = 2 s, the soil porosity almost equals the initial value, but the soil porosity increases most obviously when *T* = 10 s. Figure 15 shows the variations of *∂ϕ <sup>∂</sup><sup>t</sup>* for different wave periods. When *T* = 10 s, the values of *∂ϕ <sup>∂</sup><sup>t</sup>* at different depths reach the peak values fastest and the peak values are the biggest compared with the other wave periods. This wave period leads to the fastest wave-induced erosion in the seabed.

**Figure 13.** Distributions of the oscillatory pore-pressure and the residual pore-pressure for different wave periods: (**a**) Vertical distribution of |*Posc*|/*Pb* for different *T*; (**b**) vertical distribution of (*Pb* − |*Posc*|)/*σ* <sup>0</sup> for different *T*; (**c**) vertical distribution of *Pres* for different *T*; (**d**) vertical distribution of *Pres*/*σ* <sup>0</sup> for different *T*.

**Figure 14.** Variations of the soil porosity for different wave periods.

**Figure 15.** Variations of *∂ϕ <sup>∂</sup><sup>t</sup>* for different wave periods: (**a**) *T* = 5 s; (**b**) *T* = 10 s; (**c**) *T* = 15 s; (**d**) *T* = 20 s.

#### *6.3. Effect of Critical Concentration of Fluidized Soil Particles*

This section aims to assess the effect of critical concentration of the fluidized soil particles on the wave-induced erosion. The values of *ccr* were selected as 0.10, 0.20, 0.30, 0.40, and 0.50, respectively. Figure 16 gives the simulation results of soil porosity versus the depth for different *ccr*. It can be seen that the soil porosity increases mainly at shallow depths (within −4 m) with the growth of *ccr*. The soil at deep depths is not affected by *ccr*. The bigger the *ccr*, the more severe the soil erosion is. As shown in Figure 17, the erosion rate *∂ϕ <sup>∂</sup><sup>t</sup>* is obviously affected by *ccr* at shallow depths. Combined with Figure 17a–e, it can be seen that the peak values of *∂ϕ <sup>∂</sup><sup>t</sup>* for the selected depths increase obviously and *∂ϕ <sup>∂</sup><sup>t</sup>* reach the peak values later with the growth of *ccr*. Furthermore, the value of *∂ϕ <sup>∂</sup><sup>t</sup>* becomes negative in the later stage of the erosion process when *ccr* ≥ 0.30. It can be concluded that the bigger the *ccr*, the more remarkable the deposition effect.

**Figure 16.** Variations of the soil porosity for different critical concentrations of the fluidized soil particles.

**Figure 17.** Variations of *∂ϕ <sup>∂</sup><sup>t</sup>* for different critical concentrations of the fluidized soil particles: (**a**) *ccr* = 0.10; (**b**) *ccr* = 0.20; (**c**) *ccr* = 0.30; (**d**) *ccr* = 0.40; (**e**) *ccr* = 0.50.

#### **7. Conclusions**

In this paper, the soil erosion in silty sand seabeds induced by wave actions was numerically investigated. A kind of three-phase soil model was used in the simulation, which includes the soil skeleton, pore fluid, and fluidized soil particles. By combining the Darcy flow law, mass conservation, and mass generation equations, the wave-induced erosion process for a typical wave condition case was simulated using COMSOL Multiphysics. Then, the influences of wave height, wave period, and critical concentration of moving particles were studied. Some useful conclusions can be drawn as follows:

1. The wave-induced erosion mainly occurred at the shallow depth of the seabed. For the typical wave condition, the depth affected by the wave-induced erosion is within approximately −5.00 m. In the erosion process, the concentration of the fluidized particles increases to the critical value and then remains at a stable state within −2.00 m depth. The soil porosity and soil permeability increase significantly in the shallow seabed. The maximum values of soil porosity and soil permeability occurred at depths of about −0.50 m. It is also found that the deeper the soil, the

slower the erosion rate, and the later the peak erosion rate can reach. The numerical model proposed in this paper can be used for the analysis of the seabed coarsening phenomenon.


The seabed coarsening phenomenon commonly appears at shallow seabeds, which is because the fine particles filling in the pore space tend to move with the seepage flow under wave actions. The coarsening phenomenon of the seabed will lead to the increase of soil permeability. This is the most important effect that can significantly affect the potential and the depth of seabed liquefaction. In addition, the mechanical properties of seabed soil will also be changed with seabed coarsening. There has been no published experiment so far about the seabed erosion process induced by waves, which will be our aim in the next step.

**Author Contributions:** Methodology, Z.G.; Validation, W.Z., C.Z. and S.R.; Formal analysis, F.Y.; Writing— Original Draft preparation, C.Z. and W.Z.; Writing—Review and Editing, F.Y.; Supervision, Z.G.

**Funding:** This research was funded by National Natural Science Foundation of China (51779220, 51209183), Natural Science Foundation of Zhejiang Province (LHZ19E090003).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**




#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Simulation of a Sandy Seabed Response to Water Surface Waves Propagating on Current**

#### **Dagui Tong, Chencong Liao \*, Jinjian Chen and Qi Zhang**

State Key Laboratory of Ocean Engineering, Department of Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China; tdg123147@sjtu.edu.cn (D.T.); chenjj29@sjtu.edu.cn (J.C.); zhudoufans@sjtu.edu.cn (Q.Z.)

**\*** Correspondence: billaday@sjtu.edu.cn; Tel.: +86-135-6479-6195

Received: 25 June 2018; Accepted: 17 July 2018; Published: 20 July 2018

**Abstract:** An integrated numerical model is developed to study wave and current-induced seabed response and liquefaction in a flat seabed. The velocity-inlet wave-generating method is adopted in the present study and the finite difference method is employed to solve the Reynolds-averaged Navier-Stokes equations with *k*-ε turbulence closure. The model validation demonstrates the capacity of the present model. The parametrical study reveals that the increase of current velocity tends to elongate the wave trough and alleviate the corresponding suction force on the seabed, leading to a decrease in liquefaction depth, while the width of the liquefaction area is enlarged simultaneously. This goes against previous studies, which ignored fluid viscosity, turbulence and bed friction.

**Keywords:** wave-current-seabed interaction; RANS equations; *k*-ε model; current velocity; seabed liquefaction

#### **1. Introduction**

Water waves and currents coexist in the ocean environment, and are major loads acting on the seabed and offshore structures. There have been many reports on wave-induced structure failures [1–3]. To date, numerous studies have been carried out to explore wave–seabed interactions [4–9]. However, relevant research on wave–current–seabed interactions (WCSI) is scarce, and far from being understood. To enhance the knowledge of WCSI and make use of this knowledge in the practice of coastal and offshore engineering, studies concerning WCSI are still needed.

In the existing research on WCSI, numerical and analytical methods have been widely adopted to study wave and current-induced seabed responses. For wave–current interactions, the analytical solution of Hsu et al. [10], based on potential flow theory, is widely used, and computational fluid dynamics (CFD) methods are also employed by some researchers. Biot's theories for poro-elastic media, i.e., the quasi-static (QS), partial dynamic (*u-p*) and fully dynamic (*u-w*) theories [11], solved using the finite element method (FEM), have been used in previous studies to govern the seabed response. The inertia effects of both soil skeleton and pore fluid are excluded in QS model, and both are included in the *u-w* model. In *u-p* approximation, the inertia effect of pore fluid is ignored. Based on the results of Ulker et al. [11] and Cheng and Liu [12], Sumer [5] summarized that, for most engineering problems, both inertia effects could be neglected, particularly when involving fine sediments (silt and fine sand) or dealing with liquefaction processes. When excessive pore pressure overcomes the self-weight of seabed soil, seabed liquefaction may occur.

To the best of the authors' knowledge, Ye and Jeng [13] were the first to incorporate a third-order approximation solution concerning wave–current interactions [10] into a FEM soil model with *u-p* approximation [14]. The inertia effect of pore fluid was ignored, and the magnitude and direction of current velocity affect the seabed response significantly, especially the liquefaction depth. In recent years, a Finite Volume Method (FVM) -based numerical solution of Navier-Stokes equations has been widely used to describe fluid–solid interaction [15,16]. Instead of analytical approximation, Zhang et al. [17–19] adopted the Reynolds-averaged Navier-Stokes (RANS) equations with *k-ε* turbulence model solved by FVM to calculate the dynamic loading under wave–current interaction, and used the internal-wave-maker method [20] to generate water waves. Based on Biot's QS theory for poro-elastic mediums, Wen and Wang [21] explored the response of a two-layer seabed using the approximation of Hsu et al. [10]. Then, using the same soil model, Wen et al. [22] explored the seabed response to the combined short-crested wave and current loading. Zhang et al. [23] further extended this work to enclose the fully dynamic behavior of a seabed using Biot's fully dynamic theory [24] using the framework of Hsu et al. [10]. Recently, Yang and Ye [25] explored the residual seabed response and progressive liquefaction in a loosely-deposited seabed under wave and current loading by integrating the loading approximation of Hsu et al. [10] and a plastic soil model [26].

To overcome time and memory overconsumption of numerical simulations, analytical solutions have been proposed to explore seabed response to combined wave and current loading. Zhang et al. [27] proposed an analytical solution by integrating the third-order approximation of Hsu et al. [10] and Biot's QS theory [28]. Liu et al. [29] then extended this research to include the inertia effect of a soil skeleton with *u-p* approximation [14], which considers the acceleration of the soil skeleton. Further, Liao et al. [30] developed a new analytical approximation to study the fully-dynamic soil response, and parametrically studied the effects of wave, current, and soil characteristics on seabed response.

In summary, it can be concluded that the third-order approximation of Hsu et al. [10] concerning wave–current interactions is widely used, in both previous numerical and analytical studies. This approach utilizes the assumption of a steady uniform current and inviscid potential flow, and thus has limitations, but provides insights. As a matter of fact, a viscous water flow with potential turbulent motion and wave energy dissipates during wave propagation. Therefore, reliable simulations need to consider the effects of fluid viscosity and turbulence.

As mentioned before, Zhang et al. [17,18] developed a FVM-solved numerical model to study wave–current interactions, including turbulent motions. In their model, the wave is generated using an internal wave-maker method [20], and a corresponding "sponge layer" is used at both lateral boundaries to help eliminate wave reflections.

In the present study, a new numerical model is proposed to simulate the seabed response to combined wave and current loading. It consists of a fluid sub-model and a seabed sub-model. In the fluid sub-model, RANS equations with a *k*-*ε* turbulence enclosure are utilized to govern the wave and current-induced fluid motion with FLOW-3D v11.2 (Flow Science, Inc., Santa Fe, New Mexico, USA). Unlike the work of Zhang et al. [17–19], the wave is generated at the inlet boundary, with only one sponge layer at the outlet boundary. The finite difference method (FDM) is then used to solve the fluid motion. In the seabed sub-model, Biot's QS theory is employed to explore the seabed response using COMSOL Multiphysics 5.2 (COMSOL Inc., Burlington, MA, USA), including pore pressure, effective stresses and liquefaction, following the aforementioned summary of Sumer [5].

#### **2. Methods**

The model utilized is composed of fluid sub-model and seabed sub-model and the sub-models are integrated with the one-way coupling method (i.e., the wave pressure calculated in wave sub-model is introduced into the seabed sub-model to analyze the seabed response). The governing equations of both sub-models, the required initial and boundary conditions, and the validation of the model are described below. The WCSI is illustrated in Figure 1. A water wave train with wavelength of *Lw* (m) propagates along with an existing water current. At the outlet boundary, a sponge layer of at least one wavelength in length is set to eliminate the wave reflection. At the seafloor, a nonslip boundary can be used to simulate the fluid motion more realistically. At the air–seawater interface, the volume of fluid (VOF) method [31] is used to capture the free surface elevation.

**Figure 1.** Definition of wave–current–seabed interactions.

#### *2.1. Fluid Sub-Model*

For incompressible Newtonian fluid motion, the mass and momentum conservations are expressed with Einstein summation convention:

$$\frac{\partial \langle u\_{fi} \rangle}{\partial x\_i} = 0,\tag{1}$$

$$\frac{\partial \langle u\_{fi} \rangle}{\partial t} + \langle u\_{f\bar{j}} \rangle \frac{\partial \langle u\_{fi} \rangle}{\partial \mathbf{x}\_{\bar{j}}} = -\frac{1}{\rho\_f} \frac{\partial p\_f}{\partial \mathbf{x}\_{\bar{i}}} + \mathbf{g}\_{\bar{i}} + \frac{1}{\rho\_f} \frac{\partial}{\partial \mathbf{x}\_{\bar{j}}} \left[ \mu \left( \frac{\partial \langle u\_{fi} \rangle}{\partial \mathbf{x}\_{\bar{j}}} + \frac{\partial \langle u\_{f\bar{j}} \rangle}{\partial \mathbf{x}\_{\bar{i}}} \right) \right] - \frac{\partial \langle u'\_{f\bar{i}} u'\_{f\bar{j}} \rangle}{\partial \mathbf{x}\_{\bar{j}}},\tag{2}$$

where *uf i uf j* (*i, j* = 1, 2) and *pf* are the mean velocity (m/s) and pressure (Pa), respectively; *xi xj* is the Cartesian coordinate (*i* = 1, 2); *ρ<sup>f</sup>* is the fluid density (kg/m3); *gi* is the gravitational acceleration (m/s2) and *<sup>μ</sup>* is the molecular viscosity (Pa·s).

The turbulence influences on the mean flow field are characterized by the Reynolds stress tensor:

$$-\rho\_f \langle u'\_{fi} u'\_{f\bar{j}} \rangle = \mu\_l \left[ \frac{\partial \langle u\_{fi} \rangle}{\partial x\_j} + \frac{\partial \langle u\_{f\bar{j}} \rangle}{\partial x\_i} \right] - \frac{2}{3} \left( \rho\_f k + \mu\_l \frac{\partial \langle u\_{\bar{l}} \rangle}{\partial x\_i} \right) \delta\_{\bar{l}j\prime} \tag{3}$$

where *μ<sup>t</sup>* = *C<sup>μ</sup> k*2 *<sup>ε</sup>* is the turbulent viscosity (TKE, Pa·s), *<sup>k</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> *u f iu f i* is the turbulent kinetic energy (m2/s2), and *δij* is the Kronecker delta. The dissipation rate of TKE (*ε*) is defined as:

$$
\varepsilon = \frac{\mu}{\rho\_f} \langle \left( \frac{\partial u\_i'}{\partial x\_j} \right)^2 \rangle. \tag{4}
$$

Finally, the *k*-*ε* turbulence closure is expressed as follows:

$$\frac{\partial \mathbf{k}}{\partial t} + \langle u\_{fj} \rangle \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \left( \mu + \frac{\mu\_{l}}{\sigma\_{k}} \right) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{j}} \right] - \langle u\_{i}^{\prime} u\_{j}^{\prime} \rangle \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} - \varepsilon,\tag{5}$$

$$\frac{\partial \varepsilon}{\partial t} + \langle \mu\_{\not\slash} \rangle \frac{\partial \varepsilon}{\partial \mathbf{x}\_{\not\slash}} = \frac{\partial}{\partial \mathbf{x}\_{\not\slash}} \left[ \left( \mu + \frac{\mu\_{\not\ell}}{\sigma\_{\varepsilon}} \right) \frac{\partial \varepsilon}{\partial \mathbf{x}\_{\not\choose}} \right] + \mathbb{C}\_{1\epsilon} \frac{\varepsilon}{k} \mu\_{\not\choose} \left( \frac{\partial \mu\_{\not\choose}}{\partial \mathbf{x}\_{\not\choose}} + \frac{\partial \langle \mu\_{\not\choose} \rangle}{\partial \mathbf{x}\_{\not\choose}} \right) \frac{\partial \langle \mu\_{\not\choose} \rangle}{\partial \mathbf{x}\_{\not\choose}} - \mathbb{C}\_{2\epsilon} \frac{\varepsilon^{2}}{k} \tag{6}$$

where *σk*, *σε*, *C*2*ε*, *C*2*ε*, and *C<sup>μ</sup>* are empirical coefficients determined by experiments [32]:

$$
\sigma\_{\mathbb{R}} = 1.00, \,\sigma\_{\mathbb{C}} = 1.30, \,\,\mathbb{C}\_{1\varepsilon} = 1.44, \,\,\mathbb{C}\_{2\varepsilon} = 1.92, \,\,\mathbb{C}\_{\mathbb{H}} = 0.09. \tag{7}
$$

To diminish the influence of reflected waves from the outflow boundary, a sponge layer is set next to the outlet. In the sponge layer, the RANS equations are modified as:

$$\begin{split} \frac{\partial \rho\_f \langle u\_{fi} \rangle}{\partial t} + \frac{\partial \rho\_f \langle u\_{fi} \rangle \langle u\_{fi} \rangle}{\partial x\_{\bar{j}}} &= -\frac{\partial p\_f}{\partial x\_i} + \frac{\partial}{\partial x\_{\bar{j}}} \left[ \mu \left( \frac{\partial \langle u\_{\bar{f}i} \rangle}{\partial x\_{\bar{j}}} + \frac{\partial \langle u\_{\bar{f}i} \rangle}{\partial x\_{\bar{i}}} \right) \right] + \frac{\partial}{\partial x\_{\bar{j}}} \left( -\rho\_f \langle u\_{f\bar{i}}' u\_{f\bar{j}}' \rangle \right) + \rho\_f g\_i - \\ &\qquad \rho\_f k\_d \left( \langle u\_{fi} \rangle - \langle u\_{f\bar{i}} \rangle\_{str} \right), \end{split} \tag{8}$$

in which −*kd uf i*−*uf istr* is the artificial damping force that dissipates the wave motion, *kd* is the damping coefficient (s<sup>−</sup>1) at a given distance (*lk*, m) from the starting side of the wave-absorbing layer toward the open boundary, and *uf istr* is the background stream velocity (m/s) that is exempted from damping. The coefficient *kd* is estimated using:

$$k\_d = k\_0 + l\_k \cdot \frac{k\_1 - k\_0}{d} \,\tag{9}$$

where *k*<sup>0</sup> and *k*<sup>1</sup> (*k*<sup>1</sup> ≥ *k*0) are the values of *kd* at the starting side of the sponge layer and the open boundary, respectively. The distance *lk* is a variable measured from the starting side of the wave-absorbing layer towards the open boundary. Finally, *d* is the length of the sponge layer (m). In the present study, *k*<sup>0</sup> = 0, *k*<sup>1</sup> = 1, and *d* = 2*Lw*, where *Lw* is the incident wavelength.

#### *2.2. Seabed Sub-Model*

Biot's QS theory for a poro-elastic medium [28] is adopted to govern the seabed response. For an isotropic homogeneous sandy seabed, the conservation of mass could be expressed as:

$$
\Delta p\_s - \frac{\gamma\_w n\_s \beta\_s}{k\_s} \frac{\partial p\_s}{\partial t} + \frac{\gamma\_w}{k\_s} \frac{\partial \varepsilon\_s}{\partial t} = 0,\tag{10}
$$

in which Δ is the Laplace operator, *ps* is the pore pressure in seabed, *γ<sup>w</sup>* is the unit weight of water, *ns* is the soil porosity and *ks* is the seabed permeability. For a plane strain problem, the volume strain (*εs*) and the compressibility of pore fluid (*βs*) are, respectively, defined as follows:

$$
\varepsilon\_s = \frac{\partial u\_s}{\partial x} + \frac{\partial w\_s}{\partial z},
\tag{11}
$$

$$
\beta\_5 = \frac{1}{K\_{\rm uv}} + \frac{1 - S\_{\rm r}}{P\_{\rm uvo}},
\tag{12}
$$

where (*us*, *ws*) are soil displacements in *x-* and *z-*direction, respectively, *Kw* is the true elasticity modulus of water (taken as 2 × 109 Pa in the present study), *Pwo* is the absolute water pressure, and *Sr* is the seabed degree of saturation.

Leaving out the body forces, the equilibrium equations could be expressed as follows:

$$-G\Delta u\_s + \frac{G}{(1 - 2\nu)} \frac{\partial \varepsilon\_s}{\partial x} = -\frac{\partial p\_s}{\partial x} \,\prime \tag{13}$$

$$G\Delta w\_s + \frac{G}{(1 - 2\nu)}\frac{\partial \varepsilon\_s}{\partial z} = -\frac{\partial p\_s}{\partial z},\tag{14}$$

where *G* is the shear modulus and *ν* is the Poisson's ratio.

#### *2.3. Boundary Treatment*

In the wave sub-model, water waves are generated at the inlet boundary with linear waves, and are dissipated at the outlet boundary with the Sommerfeld radiation method [33]. Before the outlet boundary, a sponge layer of 2*Lw* long is applied to eliminate wave reflection. At the seabed surface, a no-slip boundary is applied

$$
\mu\_{fi} = \frac{\partial p\_f}{\partial z} = 0.\tag{15}
$$

in which *z* is the vertical coordinate.

The VOF method [31] is introduced to capture the free surface elevation.

In the seabed sub-model, the two lateral boundaries and the seabed bottom are set as fixed impermeable boundaries:

$$u\_s = w\_s = \frac{\partial p\_s}{\partial n} = 0,\tag{16}$$

in which *n* is the normal vector to each boundary. At the seabed surface, wave pressure (*pwv*) is applied to realize the coupling between the sub-models:

$$p\_{wv} = p\_f - \gamma\_w h\_{w\prime} \tag{17}$$

in which *γ<sup>w</sup>* is the unit weight of water and *hw* is the water depth.

#### *2.4. Numerical Scheme*

The parameters used in the present study are shown in Table 1. The incident wave is assumed to be a linear wave with wave period (*T*) of 8 s and wave height (*H*) of 3 m in 10-m-deep water. The wavelength (*Lw*) is iteratively calculated by:

$$L\_w = \frac{\lg T^2}{2\pi} \tanh(\frac{2\pi}{L\_w} h\_w). \tag{18}$$

where *g* is the gravitational force.


**Table 1.** Input data of the present study.

Wave steepness is δ = *H*/*Lw* = 0.042. To parametrically study the effect of the current velocity on seabed response, a series of current velocities from 0 to 1 m/s with a gradient of 0.25 m/s are adopted. As has been summarized by Jeng [34], the marine sediments are usually not fully saturated and have degrees of saturation very close to unity. Hence, in the present study, the seabed is considered to be unsaturated coarse sand (*Sr* = 0.985) with an isotropic permeability of 1.0 × <sup>10</sup>−<sup>4</sup> m/s. The shear

modulus (*G*), Poisson's ratio (*ν*) and porosity (*ns*) are set as 1.0 × <sup>10</sup><sup>7</sup> N/m2, 0.333 and 0.3, respectively. As an elastic seabed, the Young's modulus (*E*) is calculated by:

$$E = \mathfrak{Z}G(1+\nu). \tag{19}$$

To reach acceptable results, the model length needs to be at least two times wavelength (*Lw*) to diminish the influence of fixed boundary, as suggested by Ye and Jeng [13]. Hence, in the present model, the seabed length is set as 3*Lw* = 213 m along with a seabed thickness of 30 m. Correspondingly, the wave model length is set as 5*Lw* in which the downstream 2*Lw* long region is set as a sponge layer to minimize wave reflection.

The one-way coupling in this study is realized by introducing the wave pressure calculated from the wave sub-model at a given time to the seabed sub-model and letting the pore pressure at the seabed surface equal the wave pressure, as displayed in Equation (17). Eventually, the wave and current-induced seabed response is captured within the Biot's equations (Equations (10), (13), and (14)).

In the wave sub-model, the whole domain, including the sponge layer, is discretized into 460,850 quadrilateral cells with an element size of *H*/30 = 0.1 m. The finite difference method is used to solve the wave motion with an output data interval of *T*/40 = 0.2 s. In the seabed sub-model, the seabed consists of 159,600 quadrilateral elements with element size of 0.2 m, and is solved by FEM.

#### **3. Results**

#### *3.1. Model Validation*

In this section, the validity of the fluid sub-model will be examined against the available experimental data in the literature. The seabed sub-model has been validated using the experimental data of Tsui and Helfrich [35], Liu et al. [36], and the analytical solution of Hsu and Jeng [37], by the authors [38–40]. For example, Figure 2, modified from Tong et al. [38], displays the pore pressure response to wave loading in the experiments of Liu et al. [36] and the simulation results. It could be seen from Figure 2 that the present model reaches a good agreement with the experimental data. For more details on the model validation, readers can refer to the authors' previous work [38–40].

**Figure 2.** Validation of time series of pore pressure against experiment data of Liu et al. [36] at depths (**a**) *z* = −6.7 cm and (**b**) *z* = −26.7 cm (Adapted from Tong et al. [38]).

In this study, the experimental data of Umeyama [41] are adopted to validate the fluid sub-model. In the experiments, the wave interaction with a following current is studied with a wave flume of 25 m × 0.7 m × 1 m. The mean water depth (*hw*), current velocity (*vc*) and wave period (*T*) are kept as 30 cm, 8 cm/s, and 1 s, respectively, while three wave heights are used in the experiments. In this study, the time series of free surface elevation of wave height *H* = 3.09 cm are adopted to validate the fluid sub-model.

Figure 3 displays the time series of free surface elevation from the experiment and the present model, in which *t* is the time, *η* is the relative wave profile, i.e. the difference between free surface elevation (*z*) and still water level (*hw*), and *A* is the amplitude of the incident wave. It is seen that the present model reaches a good agreement with the experiment, in terms of wave period and amplitude, except for a slight discrepancy between the simulated and experimental results.

**Figure 3.** Validation of the fluid sub-model against the experiment data of Umeyama [41].

#### *3.2. Hydrodynamics of WCSI*

The previous studies on WCSI scarcely discussed the hydrodynamics of wave-current interaction as most of which directly use the analytical solution of Hsu et al. [10]. In this subsection, the free surface elevation and wave pressure on the seabed is shown and discussed. Among the previous studies, Zhang et al. [17] adopted the FVM-solved RANS equations with a *k-ε* turbulence closure scheme to simulate wave–current interactions with the internal wave-maker method. In the present study, based on the same governing equations, we adopted FDM to solve the problem with an incident wave generated at the inlet boundary.

Figure 4 displays the time series of free surface elevation and wave pressure on the seabed surface (location O is the midpoint on the seabed surface). It is seen that the current leads to a narrow steep crest and a flat trough of free surface elevation. As the current velocity goes up, the peak value of free surface elevation increases along with a decrease in the magnitude of the trough. Similar phenomena can also be found for the time development of wave pressure. It is known that the intensity of the wave and current-induced seabed liquefaction is dependent on the magnitude of negative wave pressure [5], i.e., negative wave pressure with larger magnitude leads to higher liquefaction potential. Therefore, it can be concluded that the existence and increase of current velocity would moderate the seabed liquefaction, which is in agreement with Zhang et al. [17]. However, this is in contrast to the result of those using the analytical solution of Hsu et al. [10] due to the assumption of uniform current omitting the effect of shear stress of the seabed. This demonstrates the necessity of considering the fluid viscosity and bed shear stress in WCSI.

**Figure 4.** Time series of (**a**) free surface elevation and (**b**) wave pressure at location O for various current velocities.

#### *3.3. Seabed Response*

Under the action of wave and current loading, excessive pore pressure will be generated in the seabed. It has been well recognized that the negative pore pressure is responsible for the wave-seabed liquefaction [4]. Correspondingly, the wave and current-induced seabed response, including the effective stresses, shear stress and pore pressure, will be presented in this subsection when the negative pore pressure reaches its maximum.

Figure 5 depicts the spatial distribution of wave and current-induced pore pressure (*ps*) and displacements (*uw*, *ws*) when magnitude of negative wave pressure at the midpoint of seabed surface reaches its maximum (*t*/*T* = 8.05) with current velocity of 0.5 m/s. As shown in the figure, seabed response to three waves is observed with an attenuation of magnitude from the inlet to the outlet. This is different from the result of Ye and Jeng [13] (Figure 6) due to the inclusion of fluid viscosity and bed friction in the present wave sub-model, which leads to the dissipation of wave energy during propagation. Besides, it is seen that the amplitude of the negative pore pressure in the horizontal plane is −12.21 kPa, which is 1.32 kPa larger than that of the positive pore pressure (10.89 kPa). Similar phenomenon is found in the distribution of vertical displacement. Ye and Jeng [13] have also shown similar phenomena in terms of pore pressure and vertical effective stress, using the analytical solution of wave–current interactions.

The effect of current velocity on seabed response has been extensively explored in the previous studies, and it is concluded that an increase of current velocity would intensify the seabed response within the analytical solution of Hsu et al. [10] on wave–current interactions. However, when the viscosity and turbulence are taken into account, Zhang et al. [18] revealed that the increase of current velocity would lead to reduction of the amplitude of pore pressure. In the present study, as shown in Figure 4, the amplitude of the positive wave pressure increases with the current velocity, while that of the negative wave pressure shows a converse trend, therefore this would lead to a different result on seabed response.

Figure 7 depicts the vertical distribution of wave and current-induced pore pressure (*ps*), effective stresses (*σ <sup>x</sup>*, *σ <sup>z</sup>*) and shear stress (*τxz*) right beneath the midpoint (Figure 5a) when the magnitude of negative pore pressure reaches its maximum, in which *h* is the thickness of the seabed. It can be observed that the magnitude of pore pressure drops down first and then increases slightly with the seabed depth, while the vertical effective stress (*σ <sup>z</sup>*) has a contrary trend. As for the horizontal effective stress (*σ <sup>x</sup>*) and shear stress, they both have a rather small magnitude throughout the seabed depth. Particularly, it can be seen that the existence and increase of current velocity lead to a decrease of the magnitude of pore pressure in the shallow seabed depth (*z*/*h* < 0.3), as well as effective stresses and shear stress. In the deep soil range (*z/h* > 0.5), it could be observed that the increase in current velocity tends to decrease the magnitude of the seabed response. Particularly, there exists a slight increase of pore pressure when seabed depth *z*/*h* increases from 0.7 to 1.0. This should be induced by

the fixed impermeable bottom boundary, which restricts the seepage and soil displacements near the seabed bottom.

**Figure 5.** Seabed response to combined wave and current loading when current velocity *vc* = 0.5 m/s.

**Figure 6.** Seabed response to combined wave and current loading when *vc* = 1 m/s in Ye and Jeng [13].

**Figure 7.** Vertical distributions of minimum wave-induced pore pressure with various current velocities.

#### *3.4. Seabed Liquefaction*

There have been several liquefaction criteria proposed in the previous studies to estimate the liquefaction potential under wave (current) loadings. In the present study, the liquefaction criterion proposed by Zen and Yamazaki [42] is adopted to evaluate the liquefaction potential in seabed:

$$-\left(\gamma\_{\mathfrak{s}} - \gamma\_{\mathfrak{w}}\right)z \le p\_{\mathfrak{s}} - p\_{\mathfrak{b}} \tag{18}$$

in which *γ<sup>s</sup>* is the unit weight of seabed soil and *pb* is the wave pressure on the seabed surface. In the previous studies, this liquefaction criterion has been widely adopted to estimate the wave-induced seabed liquefaction potential around pipelines [43,44], breakwaters [45] and pile foundations [38,39,46]. In this study, the soil unit weight is taken as 1.8*γw*.

Figure 8 illustrates the seabed liquefaction depth (*dl*) when the magnitude of negative wave and current-induced pore pressure reaches its maximum. It is seen that the maximum liquefaction depth reaches nearly 0.9 m when current velocity is zero. When there is a current, it can be seen that the liquefaction depth displays a decrease, however, with an increase in width of the liquefaction area. This corresponds to the elongation effect of current on wave trough as illustrated in Figure 4.

**Figure 8.** Maximum liquefaction depth with various current velocities.

Figure 9 presents the variations of maximum liquefaction depth (*dl*) and width (*wl*) with current velocity (*vc*). It is shown that with the increase of current velocity from 0 to 1.0 m/s, the liquefaction depth drops down. However, the width of the liquefaction area increases with the increase of current velocity.

**Figure 9.** Maximum liquefaction depth and width around location O with various current velocities.

#### **4. Discussion**

The CFD method is used in the present study to generate waves and simulate the wave–current interactions, and the RANS equations with *k-ε* turbulence model are taken as the governing equations with the VOF method to describe the free surface motion. Thus, the viscosity and turbulence of wave and current motion could be enclosed in this study in comparison with the analytical approximation of Hsu et al. [10] which simplified the fluid motion as inviscid and irrotational potential flow. Based on these governing equations, Zhang et al. [17–19] took FVM to solve the wave–current interactions with the internal wave maker method, in which sponge layers are set at both ends of the numerical flume to eliminate the wave reflection. In the present study, the FDM method is used to solve the governing equations rather than FVM. The inlet velocity method is used to generate the wave train and current with only one sponge layer at the outflow boundary. Thus, computation time and memory consumption could be saved in the present study in comparison with the model of Zhang et al. [17–19].

Based on the analytical approximation of Hsu et al. [10], the previous studies found that the increase of current velocity would aggravate the seabed liquefaction depth. However, when the RANS equations with *k-ε* turbulence model are taken to calculate the wave loading, a contrary trend is found, both in Zhang et al. [18] and the present study (Figures 6 and 7). In Hsu et al. [10], potential flow theory is adopted and the current is considered to be uniform. Hence, it goes against the fact that the current velocity near the seafloor should be quite small, even negligible due to the non-slip boundary. Correspondingly, the CFD method (RANS equations with *k-ε* turbulence model) is more reasonable when simulating wave-current interaction.

#### **5. Conclusions**

The oscillatory seabed response to combined wave and current loading is numerically explored in this paper. The FDM-solved RANS equations with *k-ε* turbulence closure and velocity-inlet wave maker are employed to simulate the wave–current interactions. Biot's QS model for a poro-elastic medium is adopted to govern the seabed response. The conclusions could be drawn as follows:

(1) The existence of current elongates the wave trough and meanwhile leads to a short wave crest. The wave energy attenuates with wave propagation, leading to magnitude attenuation of seabed response.

(2) The increase of current velocity intensifies the positive wave pressure on the seabed surface, and moderates the negative wave pressure. In the shallow seabed, the seabed response is alleviated with the current velocity while it intensifies in the deep soil range.

(3) The seabed liquefaction depth decreases with the current velocity, while the width of the liquefaction zone increases with the current velocity. This corresponds to the effect of current velocity on wave profile and wave pressure.

**Author Contributions:** Methodology, C.L.; Validation, D.T. and Q.Z.; Formal Analysis, D.T.; Writing-Original Draft Preparation, D.T.; Writing-Review and Editing, J.C.; Supervision, C.L.

**Acknowledgments:** The authors are grateful for the financial support from the National Science Foundation of China (Grant No. 41727802, No. 51678360, No. 41602282).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Effects of Principal Stress Rotation on the Fluid-Induced Soil Response in a Porous Seabed**

**Zhengxu Li 1,†, Dong-Sheng Jeng 1,\*,†, Jian-Feng Zhu 2,† and Hongyi Zhao 3,†**


Received: 22 March 2019; Accepted: 17 April 2019; Published: 28 April 2019

**Abstract:** Principal stress rotation (PSR) is an important feature for describing the stress status of marine sediments subject to cyclic loading. In this study, a one-way coupled numerical model that combines the fluid model (for wave–current interactions) and the soil model (including the effect of PSR) was established. Then, the proposed model was incorporated into the finite element analysis procedure DIANA-SWANDYNE II with PSR effects incorporated and further validated by the experimental data available in the literature. Finally, the impact of PSR on the pore-water pressures and the resultant seabed liquefaction were investigated using the numerical model, and it was found that PSR had a significant influence on the seabed response to combined wave and current loading.

**Keywords:** Principal stress rotation; dynamic loading; wave (current)-induced soil response; seabed liquefaction

#### **1. Introduction**

Recently, the physical processes of fluid–seabed interactions have attracted great attention from coastal and geotechnical engineers because of the growth in human exploration and development of offshore projects. Seabed instability due to cyclic loading, such as waves, currents, and earthquakes, is one of the main concerns of offshore geotechnical engineers involved in the design of offshore infrastructure.

It has been well known that dynamic wave pressure generated by natural hydrodynamic loading on the sea floor further induces pore-water pressure and stresses in the seabed. When the pore pressure accumulates and reaches a certain level, the effective stresses vanish and lead to soil instability as a consequence of the movement of soil particles [1,2]. Therefore, an accurate prediction of the soil response, including pore-water pressure, effective stresses, and shear stresses, is important for the design of offshore infrastructure.

On the basis of laboratory and field measurements, two mechanisms for the wave-induced soil response have been developed and reported in the literature [3–5], namely oscillatory and residual mechanisms. The first mechanism is the result of oscillatory excess pore-water pressures and accompanied by the attenuated amplitude and phase lag of pore-water pressure changes [6]. The second mechanism is the build-up of excess pore pressures caused by the contraction of cyclic loading [7]. As reported in Jeng and Seymour [8], the oscillatory mechanism dominates the process of

liquefaction in the case of a longer wave period and small amplitude, while the residual mechanism dominates the whole process for a wave with a short wave period or large wave amplitude.

Numerous studies for the wave-induced oscillatory soil response have been carried out since the 1970s. For example, on the basis of Biot's consolidation theory [9], Yamamoto et al. [6] derived a closed-form analytical solution for the wave-induced oscillatory soil response in an infinite seabed. The scope of this framework has been further extended to a seabed of finite thickness, a layered seabed, an inhomogeneous seabed with variable permeability, and a cross-anisotropic seabed [5]. Later, several analytical solutions were proposed that incorporated dynamic soil behavior; these models include the partial dynamic (*u* − *p*) [10] and full dynamic (FD) models [11,12]. Jeng and Cha [11] investigated the applicable range of different approximations with two non-dimensional parameters and soil permeability. This applicable range was reexamined by Ulker and Rahman [12] for different soils.

In addition to analytical approximations, several numerical models for the wave-induced oscillatory soil response for more complicated cases have been developed and applied to different offshore infrastructures. For example, Jeng and Lin [13] established a finite element model (FEM) that considers the effects of variable permeability and the shear modulus. Later, FEM models were developed for cases with different offshore infrastructures, for example, breakwaters [14], pipelines [15], and mono-pile foundations [16].

In the literature, numerous investigations of the wave-induced residual soil response are available. Using the results of direct shear tests [17], Seed and Rahman [7] proposed a 1D approximation with a source term for pore-water pressure generation. Following this framework, several analytical solutions and numerical solutions for wave-induced residual liquefaction were proposed [4,5,8]. Recently, Jeng and Zhao [18] proposed a new definition of the source term by considering the instant oscillatory shear stress; then, they extended the 1D model to two dimensions and applied it to the case of a submarine pipeline [19]. The aforementioned works were based on an inelastic model with a source term. Adopting the model proposed by Sassa et al. [20] and including the dissipation of pore-water pressures in the source term, Liao et al. [21] extended the model to two-dimensional cases. In addition to the inelastic models with a source term, a poro-elastoplastic model (DIANA-SWANDYNE II) was established by Chan [22] for earthquake-induced liquefaction by adopting the Pastor–Zienkiewicz Mark-III (PZIII) model [23]. This model was modified and applied to the problem of wave–seabed interactions around marine infrastructures, such as pipelines and breakwaters [24,25].

In natural ocean environments, the co-existence of waves and currents is a common physical phenomenon, and their interaction is an important topic in coastal and ocean engineering. The presence of a current in propagating waves directly changes the flow field and causes further changes to the soil response. On the basis of the analytical solution for wave–current interactions [26], Ye and Jeng [27] were the first to investigate the wave (current)-induced oscillatory soil response in a porous seabed. Following a similar framework, Wen et al. [28] further considered the case of a submarine pipeline by using the commercial software ABAQUS. Several analytical solutions based on different soil behaviors, such as quasi-static, partial dynamic, and full dynamic models, have been developed to describe the soil response to combined wave and current loading [29]. All of these approaches are based on the third-order analytical approximation for wave–current interactions [26]. Using the numerical model for wave–current interactions [30], Zhang et al. [31] further investigated the wave (current)-induced oscillatory pore pressures in a porous seabed. Recently, Liao et al. [32] proposed coupling models for residual seabed liquefaction subject to combined wave and current loading. Although numerous theoretical studies have been carried out since 2012, only two experimental studies for the wave (current)-induced oscillatory soil response are available in the literature [33,34]. These studies were used for the validation of the present model.

None of the aforementioned investigations have considered the effects of principal stress rotation (PSR) in a marine seabed, although the continuous rotation of principal stresses is an essential feature of soil's dynamic response to cyclic loading. Unfortunately, because pure PSR is assumed, this process cannot be captured by a conventional elastoplastic model without changing the cyclic deviatoric stress amplitude of the plastic strain [35]. Several experimental results have confirmed that plastic strains are generated merely by altering the principal stress orientation in both monotonic and cyclic rotational shear tests [36,37]. On the basis of the generalized plasticity theory Zienkiewicz and Morz [38], as the first attempt, Sassa and Sekiguchi [35] developed a modified version of PZIII model by considering the effects of principal stress orientation. Their model defines a new major principal stress angle parameter (Φ) that replaces potential plastic functions, the loading functions, and the plastic modulus. However, as Zhu et al. [39] pointed out, Sassa's model [35] also has deficiencies; for example, it does not account for out-of-plane stress, which is a critical parameter in the determination of plastic flow conditions. Furthermore, the reloading effect is not considered in their model. Recently, Zhu et al. [39] proposed a modified constitutive model in which both the PSR and the out-of-plane stress are taken into consideration within the generalized plasticity theory framework. In contrast to Sassa's model [35], this model was built to consider previous events during the reload by adding a discrete memory factor, and stress invariants were added at the same time to complete the optimization of the model. However, Zhu et al. [39] only considered linear wave theory in their model. The effects of PSR on the soil response to combined wave and current loading have not been reported in the literature.

In this paper, the constitutive model proposed by [39] is adopted, and the impact of PSR is included to examine the wave (current)-induced soil response in a sandy seabed. The theoretical model for both the flow model (wave–seabed interactions) and the seabed model (with PSR) are outlined first. The validation of the present model by both wave flume tests [33,34] and centrifugal tests [40] is then described. Finally, the results of the parametric study are reported to examine the effects of PSR with combined wave and current loading.

#### **2. Theoretical Models**

The present model consists of two submodels: flow and seabed submodels. A one-way coupling between the two different models is employed by the pressure continuity at the interface between the fluid and seabed domains. The fluid model is used to obtain the flow characteristics, such as wave motion, velocity field, wave pressures, etc. In the present model, the continuity of pressures is used to link the two submodels; that is, the dynamic fluid pressure is extracted and interpolated on the grid points of the solid model interface and serves as the pressure boundary condition for the seabed model. This approach has been commonly used by previous researchers [6,12,14,41,42].

The present study is based on the one-way coupling approach. Although one-way coupling has been widely used in the past and may still be effective in some cases, there are more recent approaches that effectively represent the water and bottom sediment coupled dynamics [43–46]. For example, Ran et al. [43] proposed an incompressible smoothed particle hydrodynamics (ISPH) scour model for movable bed dam break flows. Wang et al. [44] presented an ISPH simulation of scour behind seawall due to continuous tsunami overflow. Manenti et al. [45] adopted SPH model to investigate the Vajont disaster and compare their numerical simulation with 2D experiments. Wang et al. [46] further adopted their model [44] to 3D ISPH erosion model for flow passing a vertical cylinder. The technique could be further adopted to the present problem in the future.

The problem considered in this study is depicted in Figure 1. In the computation domain, the seabed thickness is *h*, and the water depth is *d*. The ocean wave propagates in the *x*-direction, while the *z*-axis is oriented upward from the seabed surface. The direction of the current can be the same as or opposite of the direction of wave propagation.

**Figure 1.** Sketch of wave (current)–seabed interaction.

#### *2.1. Flow Model*

Recently, the open-source code OpenFOAM has become widely used for the simulation of various coastal/ocean engineering problems; for example, Waves2FOAM and IHFOAM have been used to study wave generation [47,48], wave–structure interactions, and other coastal engineering processes [49]. In this study, IHFOAM was adopted to describe the wave–current interactions. Basically, IHFOAM solves three-dimensional Volume-Averaged Reynolds-Averaged Navier–Stokes (VARANS) equations for two incompressible phases (water and air) using a finite volume discretization and volume of fluid (VOF) method. The governing equations, including mass conservation and momentum conservation equations, can be expressed as

$$\frac{\partial \langle u\_{fi} \rangle}{\partial x\_i} = 0,\tag{1}$$

$$\frac{\partial \rho \langle u\_{fi} \rangle}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \frac{1}{n} \rho \langle u\_{fi} \rangle \langle u\_{fi} \rangle \right] = -n \frac{\partial \langle p^\* \rangle^f}{\partial \mathbf{x}\_i} + n \rho g\_i + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu\_{eff} \frac{\partial \langle u\_{fi} \rangle}{\partial \mathbf{x}\_j} \right] - [\mathbb{C}T]\_\prime \tag{2}$$

where and *<sup>f</sup>* are Darcy's volume-averaging operator and the intrinsic averaging operator, respectively; *ρ* is the density, computed by *ρ* = *αρwater* + (1 − *α*)*ρair*, in which *α* is the indicator function defined in (4); *uf i* is the velocity vector; *n* is the porosity; *p*<sup>∗</sup> is the pseudo-dynamic pressures; *gi* is the gravitational acceleration; *μeff* is the efficient dynamic viscosity, defined as *μeff* = *μ* + *ρνturb*, in which *μ* is the molecular dynamic viscosity and *νturb* is the turbulent kinetic viscosity, given by the chosen turbulence model. In this study, the *k* − turbulence model is used. The last term in (2) represents the resistance of porous media and can be expressed as

$$\mathbb{E}\left[\mathbb{C}T\right] = A\langle u\_{fi}\rangle + B\langle \langle u\_f \rangle \rangle \langle u\_{fi}\rangle + \mathbb{C}\frac{\partial \langle u\_{fi}\rangle}{\partial t},\tag{3}$$

where the factor *C* is less significant than factors *A* and *B* (refer to [50] for the values). A value of *C* = 0.34 [kg/m3] is often applied by default [51].

Each cell in the computational domain is considered a mixture of a two-phase fluid (air and water). The indicator function *α* varies from 0 (air) to 1 (water); *α* is defined as the quantity of water per unit of volume for each cell and calculated as follows:

$$\mathfrak{a} = \begin{cases} 1, & \text{water} \\ 0, & \text{air} \\ 0 < \mathfrak{a} < 1, & \text{free surface} \end{cases} \tag{4}$$

Any variation of fluid properties, such as density and viscosity, can be represented using the indicator function *α* considering the mixture properties:

$$
\Phi = \mathfrak{a}\Phi\_{water} + (1-\mathfrak{a})\Phi\_{air} \tag{5}
$$

where Φ*water* and Φ*air* are water and air properties, respectively, such as the density of the fluid.

The fluid's movement can be tracked by solving the following advection equation [52]:

$$\frac{\partial \alpha}{\partial t} + \frac{1}{n} \frac{\partial \langle u\_{fi} \rangle \alpha}{\partial x\_i} + \frac{1}{n} \frac{\partial \langle u\_{fi} \rangle \alpha (1 - \alpha)}{\partial x\_i} = 0,\tag{6}$$

where <sup>|</sup>**u***f c*<sup>|</sup> <sup>=</sup> *min cα*|**u***<sup>f</sup>* |, *max*(|**u***<sup>f</sup>* |) , in which the default value of *cα* is 1.

The wave generation and active wave absorption in the fluid domain were implemented within IHFOAM. Several boundary conditions were introduced: (i) the inlet boundary condition allows for generating a wave according to different wave theories as well as adding different steady current flows; (ii) the outlet boundary condition applies an active wave absorption theory to prevent the re-reflection of an incoming wave; (iii) the slip boundary condition (zero-gradient) is applied on the bottom of the fluid domain and the lateral boundary of the numerical wave flume; (iv) the top boundary condition is set as the atmospheric pressure. For the details of IHFOAM, the readers can refer to Higuera et al. [48].

#### *2.2. Seabed Model*

In the literature, three different models of fluid–seabed interactions have been established on the basis of different soil behaviors: quasi-static (QS, i.e., the conventional consolidation model), partial dynamic (*u* − *p*), and full dynamic (FD) models. All three are based on Biot's porous theory [9,53]. Zienkiewicz et al. [54] proposed the *u* − *p* approximation and examined the applicable range between *u* − *p* and QS for earthquake loading. The framework has been further extended to the problem of wave–seabed interactions [10–12]. The applicable range of different models has been clarified for various soil types [11,12].

This paper establishes a two-dimensional model that considers the rotation of the principal stress axis to analyze the seabed response to combined dynamic loading due to wave–current interactions. The dynamic Biot's equation proposed by Zienkiewicz et al. [54], the *u* − *p* approximation, was adopted.

$$\frac{\partial \sigma\_x'}{\partial \mathbf{x}} + \frac{\partial \tau\_{xz}}{\partial z} = -\frac{\partial p\_\varepsilon}{\partial \mathbf{x}} + \rho \frac{\partial^2 u\_\delta}{\partial t^2} \,, \tag{7}$$

$$\frac{\partial \tau\_{xz}}{\partial x} + \frac{\partial \sigma\_z'}{\partial z} + \rho b\_\% = -\frac{\partial p\_\ell}{\partial z} + \rho \frac{\partial^2 w\_s}{\partial t^2},\tag{8}$$

$$K\_s \nabla^2 p\_\varepsilon - \gamma\_w n\_s \beta\_s \frac{\partial p\_\varepsilon}{\partial t} + K\_s \rho\_f \frac{\partial^2}{\partial t^2} \left(\frac{\partial u\_s}{\partial x} + \frac{\partial w\_s}{\partial z}\right) = \gamma\_w \frac{\partial \varepsilon\_v}{\partial t},\tag{9}$$

where *σ <sup>x</sup>* and *σ <sup>z</sup>* are the effective normal stresses in the *x*- and *z*-directions, respectively; *τxz* is the shear stress; *pe* is the pore-water pressure; *us* and *ws* represent the soil displacement in the *x*- and

*<sup>z</sup>*-directions; *bg* is the gravitational acceleration; *Ks* is the soil permeability; ∇<sup>2</sup> is the Laplace operator; *n* is soil porosity; *ρ* is the average density of a porous seabed and defined by *ρ* = *ρ<sup>f</sup> n* + *ρs*(1 − *n*), in which *ρ<sup>f</sup>* is the fluid density while *ρ<sup>s</sup>* is the solid density.

In (9), the compressibility of the pore fluid *β<sup>s</sup>* is defined as [55]

$$\beta\_s = \frac{1}{K\_w} + \frac{1 - \mathcal{S}\_r}{P\_{w0}},\tag{10}$$

where *Kw* is the true bulk modulus of the elasticity of water (which can be taken as 1.95×109 N/m<sup>2</sup> [6]), *Sr* is the degree of saturation, and *Pw*<sup>0</sup> is the absolute water pressure. When the soil is fully saturated (i.e., it is completely air-free), then *β<sup>s</sup>* = 1/*Kw* since *Sr* = 1.

The anisotropic elastic constitutive model cannot account for the directional effect of the principal stress or the dilatancy of sand. Compared with the elastic constitutive model, plastic constitutive models can more realistically simulate the stress–strain relationship of soil under dynamic load conditions and the accumulation of pore-water pressure. Therefore, Zhu et al. [39] proposed a plastic constitutive model, which was implemented in the DIANA-SWANDYNE II [22] finite element code. This code is used to analyze the seabed response of the principal stress axis to waves and ocean currents. In Zhu's plastic constitutive model [39], the loading direction vector *nij* = *∂ f* /*∂σij*, the plastic flow direction vector *mij* = *∂g*/*∂σij*, and plastic modulus *HL*(*p* , *q*, *θ*, *ψ*) are defined. For the theory of generalized plasticity, the loading function *f*(*p* , *q*, *θ*, *ψ*) and plastic potential function *g*(*p* , *q*, *θ*, *ψ*) do not need to be explicitly defined.

The elastic–plastic constitutive matrix can be expressed in tensorial notation as

$$d\sigma'\_{ij} = (D^\varepsilon\_{ijkl} - \frac{D^\varepsilon\_{ijmn} m\_{mm} n\_{stl} D^\varepsilon\_{stkl}}{H\_{L,kl} + n\_{st} D^\varepsilon\_{stkl} m\_{kl}}) d\_{\text{ckl}\prime} \tag{11}$$

where *D<sup>e</sup> ijkl* is the elastic stiffness tensor.

The loading direction vector *nij* and the plastic flow direction vector *mij* can be defined as

$$n\_{ij} = \frac{\frac{\partial f}{\partial \sigma'\_{ij}}}{\left\| \frac{\partial f}{\partial \sigma'} \right\|} + \chi \frac{\partial \psi}{\partial \sigma'\_{ij}} = \frac{\frac{\partial f}{\partial p'} \frac{\partial p'}{\partial \sigma'\_{ij}} + \frac{\partial f}{\partial q} \frac{\partial q}{\partial \sigma'\_{ij}} + \frac{\partial f}{\partial \theta} \frac{\partial \theta}{\partial \sigma'\_{ij}}}{\left\| \frac{\partial f}{\partial p'} \frac{\partial p'}{\partial \sigma'} + \frac{\partial f}{\partial q} \frac{\partial q}{\partial \sigma'} + \frac{\partial f}{\partial \theta} \frac{\partial \theta}{\partial \sigma'} \right\|} + \chi \frac{\partial \psi}{\partial \sigma'\_{ij}},\tag{12}$$

$$m\_{ij} = \frac{\frac{\partial g}{\partial \sigma'\_{ij}}}{\left\| \frac{\partial g}{\partial \sigma'} \right\|} = \frac{\frac{\partial g}{\partial p'} \frac{\partial p'}{\partial \sigma'\_{ij}} + \frac{\partial g}{\partial q} \frac{\partial q}{\partial \sigma'\_{ij}} + \frac{\partial g}{\partial \theta} \frac{\partial \theta}{\partial \sigma'\_{ij}}}{\left\| \frac{\partial g}{\partial p'} \frac{\partial p'}{\partial \sigma'} + \frac{\partial g}{\partial q} \frac{\partial q}{\partial \sigma'} + \frac{\partial g}{\partial \theta} \frac{\partial \theta}{\partial \sigma'} \right\|}\,\tag{13}$$

and

$$M\_f\left(\psi\right) = M\_{f0} - \mathcal{U}\left(\psi\right)aM\_{f0}, \qquad M\_{\mathcal{g}}\left(\psi\right) = M\_{\mathcal{g}0} - \mathcal{U}\left(\psi\right)aM\_{\mathcal{g}0} \tag{14}$$

$$M\_{f0} = \frac{18M\_{f\varepsilon}}{18 + 3\left(1 - \sin 3\theta\right)}, \qquad M\_{\%} = \frac{18M\_{\%}}{18 + 3\left(1 - \sin 3\theta\right)}\tag{15}$$

$$
\mathfrak{a}\left(\psi\right) = \mathfrak{a}\_0 + c\mathcal{U}\left(\psi\right) \tag{16}
$$

$$\mathcal{U}\left(\psi\right) = \begin{cases} 1 - \cos\left(2\psi\right) & 0 \le |\psi| \le \frac{\pi}{4} \\ 1 - \cos\left(2\left|\psi\right| - \pi\right) & \frac{\pi}{4} \le |\psi| \le \frac{\pi}{2} \end{cases} \tag{17}$$

where *χ* is he control parameter to account for the effect of PSR; *ψ* is the major principal stress angle; *α*0, a, and c are the principal stress orientation model parameters; and *Mf c* and *Mgc* are model parameters related to the stress ratio.

In addition, taking into account the effects of PSR, the plastic modulus *HL* is given as

*J. Mar. Sci. Eng.* **2019**, *7*, 123

$$H\_L = H\_0 p' \left[ 1 - \frac{\eta(\psi)}{\eta\_f^\*} \right]^4 \left[ 1 - \frac{q/p'}{M\_\circ(\psi)} + \beta\_0 \beta\_1 exp(-\beta\_0 \xi) \right] H\_{DM} \tag{18}$$

where *H*0, *β*0, and *β*<sup>1</sup> are model parameters, and

$$\eta(\psi) = \frac{q}{p'} + [1 - \mathcal{U}(\psi)] a \mathcal{M}\_{\mathbb{S}^{0, \prime}} \tag{19}$$

$$
\eta\_f^\* = (M\_{f0} - aM\_{\%}) \left( 1 + \frac{1}{\kappa\_0 + \varepsilon} \right),
\tag{20}
$$

$$
\xi^\mathfrak{x} = \int d\xi^\mathfrak{x} = \int |d\xi^\mathfrak{x}|.\tag{21}
$$

where *ξ* is the accumulated deviatoric plastic strain.

To consider the history of loading events throughout the reloading process, the discrete memory factor *HDM* is introduced in the following:

$$H\_{DM} = \left(\frac{\zeta\_{max}}{\zeta}\right)^{\gamma d} \tag{22}$$

$$\mathcal{L} = p' \left\{ 1 - \left[ \frac{1 + a(\psi)}{a(\psi)} \right] \frac{q / p'}{M\_{\mathcal{S}}(\psi)} \right\}^{1/a(\psi)} \tag{23}$$

and *γ<sup>d</sup>* is the coefficient for the discrete memory factor.

The plastic modulus *HU* for unloading is

$$H\_{lI} = \begin{cases} H\_{lI0} \left( \frac{M\_\varepsilon(\varphi)}{q/p'} \right)^{\gamma\_u}, & for \quad \left| \frac{M\_\varepsilon(\varphi)}{q/p'} \right| > 1\\ H\_{lI0}, & for \quad \left| \frac{M\_\varepsilon(\varphi)}{q/p'} \right| \le 1 \end{cases} \tag{24}$$

where *HU*<sup>0</sup> and *γ<sup>u</sup>* are original model parameters.

#### **3. Model Verification**

To validate the present model, two comparisons with previous experimental data are presented here. First, we compare the present model with the hollow cylinder apparatus (HCA) element tests [56] for the present constitutive model. Second, we compare the present model with the wave flume test for wave (current)-induced oscillatory pore-water pressures [33,34]. Third, we compare the present model with centrifugal tests [40] for the development of pore-water pressure build-up.

#### *3.1. Comparison with Hollow Cylinder Apparatus (HCA) Element Tests*

In the validation of the constitutive model with regard to PSR, the HCA elementary test is urgent in need due to its particular ability in simulation of the PSR through altering its control parameters such as axial load, torque, inner cell pressure, and outer cell pressure. Towhata and Ishihara [56] such a test with pure rotation of principal stress axis, in which the major principal stress orientation angle ranged from −*π*/4 to *π*/4 with a constant deviatoric stress of 76.7 kPa. The main model parameters are shown in Table 1. The predicted results of the present model with PSR are plotted in Figure 2, in which the test data of Towhata and Ishihara [56] are also illustrated for comparison. It can be clearly seen that the adopted constitutive model present behaves well in capturing the effect of PSR on the development of volumetric strains with the constant amplitude of deviatoric stress. Moreover, the feasibility of the constitutive model with PSR is well validated by the good agreement between the predicted results and the measured data as shown in Figure 2.

**Table 1.** Present constitutive model's parameters for comparison with the HCA (Hollow Cylinder Apparatus) element test of sand [56].


(**b**) Volume change

**Figure 2.** Comparison of the present model with the HCA tests [56] under continuous rotation of the major principal stress axis.

#### *3.2. Comparison with Laboratory Experiments for the Seabed Response to Waves and Currents*

The second verification of the present model compares the present model's results with experimental results. Qi and Gao [33] conducted a series of experiments to investigate the seabed response to different wave and current velocities around a single pile. In their experiments, a wave flume of 52 m × 1 m with a depth of 1.5 m was set over sandy soil. The depth of the sandy tank was 0.5 m. However, the data for the wave and opposite current cannot be adopted for the comparison with the numerical model because of the effect of the single pile, which was set in the middle of the experiments. Thus, only the cases of the wave alone and the wave with a following current (when waves and currents have the same direction of propagation) were used to verify the model. In their experiments, the pore-pressure build-up (i.e., residual mechanism) was not observed, i.e., the tests are in the range applicable to the oscillatory mechanism, for which the elastic model is used for comparison here.

Using the same conditions as those of the wave flume tests [34], the following values are assumed: the water depth is 0.5 m; the soil model is set below the wave flume with a length of 2.4 m; and the thickness of the saturated soil layer is 0.5 m. The wave profile at the free surface and the corresponding pore-water pressure beneath 0.1 m of the seabed surface are compared. The numerical results of the wave profile and pore pressures have an overall agreement with the experimental data, as shown in Figures 3 and 4.

In the above comparisons, the top subfigures show the water surface elevation and the bottom subfigures show the pore-water pressures at 10 cm beneath the seabed surface. As shown in Figures 3 and 4, the present model predicts the wave profiles (*η*) well, but the pore-water pressures present some differences between the numerical results and experimental data, although the general trends are in agreement.

**Figure 3.** Validation of the present model with the experimental data [33] (wave only): (**a**) water surface elevation (*η*) and (**b**) the wave-induced pore pressure in seabed (*ps*). Input data: wave height (*H*)=5 cm, wave period (*T*) = 1 s, water depth (*d*) = 50 cm, seabed thickness (*h*) = 50 cm, degree of saturation (*Sr*) = 1, Shear modulus (*G*) = 107 N/m2, Poisson's ratio (*μ*) = 0.3, soil permeability (*Ks*) = 1.88 <sup>×</sup> <sup>10</sup>−<sup>4</sup> m/s, and soil porosity (*ns*) = 0.771.

**Figure 4.** Validation of the present model with experimental data resulting from combining the wave and the following current loading [33]: (**a**) water surface elevation (*η*) and (**b**) the wave-induced pore pressure in seabed (*ps*). Input data: *H* = 5 cm, *T* = 1 s, *d* = 50 cm, *h* = 50 cm, *Sr* = 1, *G* = 107 N/m2, *<sup>μ</sup>* = 0.3, *Ks* = 1.88 <sup>×</sup>10−<sup>4</sup> m/s, *ns* = 0.771, current velocity (*U*) = 0.05 m/s.

#### *3.3. Comparison with Centrifuge Tests and Previous Numerical Model for the Seabed Response To Waves*

Sassa and Sekiguchi [40] carried out a series of geotechnical centrifuge tests to investigate the process of wave-induced liquefaction in a sandy seabed. To verify the present model, we reproduced the experimental conditions and compared the results with the centrifugal experimental data [40]. Their wave tests were all performed with a centrifugal acceleration of 50 *g* (where *g* is the gravitational acceleration). The soil bed was 200 mm in width and 100 mm in depth. The submerged unit weight of soil, *γ* , was equal to 425 kN/m3, and the wave number, *κ*(= 2*π*/*L*0), was 12.2 m<sup>−</sup>1. The corresponding wave loading intensity *p*<sup>0</sup> and cyclic stress ratio (*χ*<sup>0</sup> = *κp*0/*γ* ) were 5.0 kPa and 0.14, respectively. Other input data are listed in Table 2. In the numerical model, we converted the problem back to an environment with a gravitational acceleration of 1. As shown in Figure 5, in general, the present model has an overall agreement with the centrifuge tests [40]. By examining the comparisons closely, we observe that the present model can capture the magnitude of the maximum pore-water pressures and the time it takes to reach the maximum pore water pressures. However, there are some differences between the numerical prediction and the centrifugal tests that occur during the pore-water pressure build-up. This implies that the present model requires further improvement. However, the magnitude of the maximum pore pressures directly affects the liquefaction depth, which is more important for

practical engineering design. Therefore, the present model can provide sufficient information for engineering design.

**Figure 5.** Comparison of the distribution of excessive pore pressure between the present model (solid line) and the centrifuge tests (dashed line) [40].

In addition to the comparison with the centrifugal tests [40], the present model is also compared with Sassa's numerical model [35] in Figure 6. The phenomenon of pore pressure build-up can be observed in the first several wave cycles. After a certain wave period, *presidual* reaches its peak value and stabilizes because liquefaction occurs. In the present model, the calculated time it takes for the residual pore pressure to reach its peak is 1100 μs, which is less than that predicted by Sassa and Sekiguchi [35] (1200 μs). This is because Sassa and Sekiguchi [35] did not consider the effect of out-of-plane stress, which is important for determining the plastic flow direction [57–59].


**Table 2.** Parameters used for comparison between the centrifugal test and numerical model.
