**1. Introduction**

In recent years, to meet the continual improvement requirements in coastal transportation, the immersed tunnel has become one of the choices to fundamentally transform the transport in the region of oceans and rivers, replacing the conventional methods such as the ferry. The immersed tunnel has a history of about 100 years and it shows a good performance in reliability and applicability under complex natural dynamic loading; for example, the longest immersed tunnel in the world, the Hongkong–Zhuhai–Macao Bridge immersed tunnel. As an alternative to a bridge, the immersed tunnel has advantages in less environmental effects and no obstruction of navigation channels. Compared to bridges, the immersed tunnel is commonly constructed in a soft and loose seabed. Thus, the stability of the seabed soil surrounding the immersed tunnel in complex marine environments becomes one of the main concerns implicated in tunnel design and maintenance.

It has been recognized that the pore water pressures and stresses in seabeds are affected by the water pressures generated by the natural dynamic loading. If the pore water pressure reaches the limit value, the liquefaction could occur with the effective stress in seabed vanishing. To avoid seabed instability around the immersed tunnel, the study of seabed dynamic behaviour is necessary under the real hydrodynamic loading. Two mechanisms of wave-induced liquefaction has been figured out based on a mass of laboratory tests and field exploration [1–3], which are transient liquefaction and residual liquefaction. The transient liquefaction is motivated by the oscillatory excess pore water pressures under wave pressure vibration which usually happens with amplitude reduction and phase lag of pore pressure in seabed soil [4]. While the residual liquefaction is on the consequence of the excess pore water pressure build-up under cyclic wave loading [5]. Later, Jeng and Seymour [6] proposed a simplified approximation to predict the liquefaction process in a large seabed, and concluded that the residual mechanism is more essential under large waves while the transient mechanism is dominate in a seabed under small wave loading.

In the past 40 years, various analytical formulas have been developed and verified in regard to the seabed dynamic response [4,7,8]. However, the seabed dynamic response around the structure is difficult to be described by analytical methods. Thus, several numerical model was developed to simulate the soil behaviour around the offshore structures, such as breakwaters [9], and buried pipes [10]. However, the research of a wave-induced response around the undersea immersed tunnel is quite limited. For instance, a real case simulation of the Busan–Geoje fixed link in South Korea was conducted by Kasper et al. [11] under large waves (wave height up to 9.2 m) generated by typhoons. Nevertheless, this study did not consider the impact of the seabed liquefaction around the tunnel. Recently, [12,13] simulated the seabed transient and residual response around the immersed tunnel under wave loading based on Biot's consolidation equations neglecting the inertial terms for soil skeleton and fluid phase.

In natural ocean environments, current is another crucial component besides wave. For instance, a long-term mentoring data of the Lingding Bay, in which the Hong Kong–Zhuhai–Macau bridge tunnel located, shows the current co-exists with wave varying from bottom to surface as the consequence of the irregular semi-diurnal tide. The maximum velocity of the surface current reaches up to 1.5 m/s [14]. The interaction between current and wave is found to be able to affect the hydrodynamic properties directly and further impact the porous seabed dynamics. Ye and Jeng [15] investigated the transient response of a porous seabed under wave combined current loading firstly. Later, the current effects in the vicinity of a submarine pipeline were examined by Wen et al. [16] based on the commercial software ABAQUS. Lately, Liao et al. [17] simulated the residual seabed liquefaction under the flow field that wave and current generated simultaneously. To date, how current affect the liquefaction of seabed soil surrounding the immersed tunnel has not been examined yet to the author's knowledge.

The aforementioned numerical models mainly developed adopting the conventional methods such as the finite element method and finite difference method. In recent years, a new numerical technique has come up which uses a set of nodes instead of the conservative meshes to approximate the solution. In consequence of discretization of the partial differential equations directly on nodes instead of meshes, the meshfree method is more qualified in dealing with mesh entanglement problems and constructing the approximations with arbitrary order of continuity than the traditional mash-based numerical methods. The start point of the meshfree method was using the moving least-square (MLS) method to establish shape functions for a set of scatter nodes by Nayroles et al. [18]. After that, Belytschko et al. [19] used the Galerkin weak form to improve the diffuse element method, which formed a new element-free Galerkin method. The element-free Galerkin method has been widely used in soil mechanics problems. In addition, an innovative interpolation scheme to overcome the disadvantage of the MLS was proposed by Liu and Gu [20], which is called the point interpolation method (PIM). The original PIM method picks up the polynomial basis as the basis function which performs well for one dimensional problems. However, it is hard to determine the ranks when this method extended to multi-dimensional range, which is the result of basis function selection. In order to figure this problem, Wang et al. [21] proposed a point interpolation method based on the radial basis functions (RBF). This method maps the multi-dimensional space into one-dimensional space though a radial function, which makes choosing basis function easier. This method has been widely used

on the geomechanics problems. For example, Wang and Liu [22] simulated the Biot's consolidation process by radial PIM method. The displacement and the pore water pressure were discretised by the same shape function in space, while the fully implicit integration scheme was adopted for time discretization to avoid the spurious ripple effect. In the literatures, two different type of the RBFs are commonly used, which are the multiquadrics (MQ) [23], and the Gaussian radial basis functions [24]. In this study, the MQ was used. These functions were first under intensive research in multivariate data interpolation and used to solve the partial differential equations by Kansa [25]. Thus, this method is usually called Kansa's method, and is known as the global RBF collocation method (GRBFCM). This method has been applied to computational fluid dynamics problems such as solutions of the Navier–Stokes equation [26], natural convections in porous media [27], numerical wave tanks [28], and solid–liquid phase change problems [29]. The GRBFCM has the obvious advantages in dealing with arbitrary and complex domain, which applied in many field, such as surface fitting, turbulence analysis, neural networks, meteorology. However, ill-conditions can occur when the resolution is high. For the purpose of overcoming the difficulties of ill-conditions as mentioned above, a localization procedure to transform the dense matrices into sparse matrices was proposed. On the basis of the multiquadric RBF [30], Lee et al. [31] proposed the local RBF collocation method (LRBFCM) for the first time. This method has been applied in various fields, such as the solutions of diffusion problem [32], Darcy flow in porous media [33], water wave scattering [34], macro-segregation phenomena [35] and so on. Recently, Wang et al. [36] adopted the LRBFCM based on the Biot's consolidation equation to simulate a wave-induced seabed response around a pipeline.

In this paper, an integrated numerical model is proposed to simulate the sandy seabed dynamic response and transient liquefaction in the vicinity of an immersed tunnel under natural complex loading. The flow model is developed based on IHFOAM [37], while the seabed model is established adopting the LRBFCM with Biot's "*u* − *p*" approximation which considered the inertial term of soil skeleton. The verification of the integrated model is carried out by comparisons with the analytical solution [7] and published experimental data [38]. In this study, the effects of the current on the seabed behaviour around the immersed tunnel are examined. Parametric studies are conducted in regard of the different wave characteristics and soil properties.
