*2.4. Meshfree Model for the Seabed Domain*

The LRBFCM is adopted to solve the governing equations listed above. First, an approximation function Φ(*yi*) which can either stand for displacement or pore pressure in the "*u* − *p* 00 formulation is considered in the computational geometry. This function is composed of an arbitrarily distributed points series *yj*(*j* = 1, 2, · · · , *n*) located both in the computational domain and on its boundary. Approximate Φ = Φ(*x*) around *y<sup>i</sup>* by the RBF *χ*(*rm*) to construct a linear equation for each node *y<sup>n</sup>* as,

$$\Phi(\boldsymbol{x}) \approx \sum\_{m=1}^{K} \alpha\_m \chi(r\_m),\tag{21}$$

where *α<sup>m</sup>* is the undetermined coefficient for *χ*(*rm*), and *χ*(*rm*) is the MQ defined by,

$$
\chi(r\_m) = \sqrt{r\_m^2 + c^2} \,\,\,\,\,\,\tag{22}
$$

with *rm*, *c*, and *x<sup>m</sup>* being the Euclidean distance from *x* to *xm*, the shape parameter [30], and the positions of the *K* nearest neighbour nodes around the prescribed centre *x*<sup>1</sup> = *yn*, respectively. An algorithm based on the kd-tree is adopted to search the *K* nearest neighbour nodes [47].

Then, (21) is collocated on the *K* nearest neighbour nodes, which can be expressed as:

$$[\Phi]\_{\mathbf{K}\times\mathbf{1}} = [\chi]\_{\mathbf{K}\times\mathbf{K}} [\mathfrak{a}]\_{\mathbf{K}\times\mathbf{1}} \tag{23}$$

$$\left[\Phi\right]\_{\mathbf{K}\times\mathbf{1}} = \left[\Phi(\mathbf{x}\_1), \Phi(\mathbf{x}\_2), \dots, \Phi(\mathbf{x}\_{\mathbf{K}})\right]^T,\tag{24}$$

$$[\chi]\_{\mathbf{K}\times\mathbf{K}} = \begin{bmatrix} \chi(\|\mathbf{\boldsymbol{x}}\_1 - \mathbf{x}\_1\|) & \chi(\|\|\mathbf{x}\_1 - \mathbf{x}\_2\|) & \cdots & \chi(\|\|\mathbf{x}\_1 - \mathbf{x}\_K\|) \\ \chi(\|\|\mathbf{x}\_2 - \mathbf{x}\_1\|) & \chi(\|\|\mathbf{x}\_2 - \mathbf{x}\_2\|) & \cdots & \chi(\|\|\mathbf{x}\_2 - \mathbf{x}\_K\|) \\ \vdots & \vdots & \ddots & \vdots \\ \chi(\|\|\mathbf{x}\_K - \mathbf{x}\_1\|) & \chi(\|\|\mathbf{x}\_K - \mathbf{x}\_2\|) & \cdots & \chi(\|\|\mathbf{x}\_K - \mathbf{x}\_K\|) \end{bmatrix},\tag{25}$$

$$[\mathbf{a}]\_{K\times1} = \begin{bmatrix} a\_1 \\ a\_2 \\ \vdots \\ a\_K \end{bmatrix},\tag{26}$$

Inversion (23) gives,

$$\left[\mathfrak{a}\right]\_{\mathbf{K}\times\mathbf{1}} = \left[\chi\right]\_{\mathbf{K}\times\mathbf{K}}{}^{-1}\left[\mathfrak{d}\right]\_{\mathbf{K}\times\mathbf{1}}\tag{27}$$

Here, we assume a linear differential operator of both governing equation and the boundary condition which represented by *L*. If collocated the *L*Φ(*x*) on *x*<sup>1</sup> = *yn*, (23) can be written as:

$$L\Phi(y\_n) = \sum\_{m=1}^{K} \alpha\_m L\chi\left(r\_m\right)|\_{\mathbf{x}=\mathbf{x}\_1} \tag{28}$$

which can be expressed as the vector form:

$$L\Phi(y\_n) = [L\chi]\_{1\times K} \, [\mathfrak{a}]\_{K\times 1} \tag{29}$$

In (29) mentioned above, the *Lχ*(*rm*) is on behalf of the results of the differential operator act on the the RBF *χ*(*rm*). Then, by combining (29) and (27), it can be obtained as:

$$L\Phi(y\_n) = [\mathbb{C}]\_{1\times K} [\phi]\_{K\times 1} \tag{30}$$

$$[\mathbb{C}]\_{1\times K} = [L\chi]\_{1\times K} \left[\chi\right]\_{K\times K}^{-1} \tag{31}$$

$$L[\boldsymbol{L}\chi]\_{1\times K} = \begin{bmatrix} L\chi(\boldsymbol{r}\_1)|\_{\mathbf{x}=\mathbf{x}\_1} & L\chi(\boldsymbol{r}\_2)|\_{\mathbf{x}=\mathbf{x}\_2} & \cdots & L\chi(\boldsymbol{r}\_K)|\_{\mathbf{x}=\mathbf{x}\_K} \end{bmatrix} \tag{32}$$

Obviously, the value of the row vector [*C*]1×*<sup>K</sup>* can be obtained if all variables *L*, *χ* and *x<sup>j</sup>* are known.

For all *y<sup>n</sup>* in computational domain, the linear equations can be obtained by the above-mentioned localization procedure on either the governing equation or boundary conditions. Next, these equations can be assembled into the system matrix, as:

$$[A]\_{N \times N} [\Phi]\_{N \times 1} = [B]\_{N \times 1} \tag{33}$$

where

$$\begin{aligned} \left[\Phi\right]\_{N\times 1} &= \left[u\_1(y\_1), u\_1(y\_2), \dots, u\_1(y\_N), u\_2(y\_1), u\_2(y\_2), \dots, u\_2(y\_N), \\ &\quad p(y\_1), p(y\_2), \dots, p(y\_N)\right]^T, \end{aligned} \tag{34}$$

which, [*φ*]*N*×<sup>1</sup> is the sought solution, while [*B*]*N*×<sup>1</sup> is a column vector contributed from the external loadings. It is noticeable that (33) is a sparse system matrix which is similar to both finite difference method as well as the finite element method. In the present numerical model, the direct solver of SuperLU [48] is adopted to solve the sparse system, (33). The procedure of the LRBFCM is finished here.

In numerical strategy, it is necessary to integrate the governing equations in the time domain if the boundary conditions or the extrernal loadings are time-dependent. In the present model, the single-step time integration method of the Newmark method [49] is adopted, which could handle each time step independently when the first- and second-order time derivatives exist at the same time. More details can be found in some previous studies [50,51].
