## 2 Random Finite Element Method
### Numerical Model  
The finite element method (FEM) was employed to evaluate the bearing capacity of the skirted foundations in homogeneous soil. In accordance with the previous research by [PERSON] et al. [16], a skirted foundation with a diameter (\(D\)) of 10 m and length (\(L\)) of 10 m was modeled in the commercial finite element software ABAQUS (version 6.14) as a perfectly rigid body. The skirt thickness (\(t\)) was set to 0.03 m. Horizontal, vertical, and moment loading (\(V\), \(H\), and \(M\)) were applied to the reference point (\(RP\)) to determine the corresponding bearing capacity, with \(RP\) located at the centre of the skirted foundation lid. Figure 2 presents the details of the skirted foundation and the positive loading and displacement. \begin{table}
\begin{tabular}{c c c c} \hline
**Case** & **COV** & \(\delta_{h}\) **(m)** & \(\delta_{v}\) **(m)** \\ \hline Ani-1 & 0.1 & 50 & 4 \\ Ani-2 & 0.2 & 50 & 4 \\ Ani-3 & 0.3 & 50 & 4 \\ Ani-4 & 0.4 & 50 & 4 \\ Ani-5 & 0.5 & 50 & 4 \\ Ani-6 & 0.3 & 30 & 4 \\ Ani-7 & 0.3 & 40 & 4 \\ Ani-8 & 0.3 & 60 & 4 \\ Ani-9 & 0.3 & 50 & 2 \\ Ani-10 & 0.3 & 50 & 6 \\ Ani-11 & 0.3 & 50 & 8 \\ \hline \end{tabular}
\end{table} 
Table 1: Summary of the parameters of soil random field. In reality, skirted foundation has a circular or square section, but it was considered that plan geometry is secondary to soil strength profile in terms of the effect on bearing capacity by [PERSON] and [PERSON] [2]. Therefore, a two-dimensional (\(2D\)) plane strain condition was assumed in this study to allow comparison with existing analytical solution [2], as well as improving calculation efficiency with satisfactory accuracy [16]. The soil domain was modeled as a rectangular region with dimensions \(7D\) in length and \(5D\) in width (\(70\times 50\) m) to minimize the effect of boundaries. Lateral boundaries were constrained in the horizontal direction, while the bottom boundary was fixed in both horizontal and vertical directions. No constraint was applied in the top boundary. The soil domain was meshed into a square grid with dimensions of \(1\times 1\) m; however, in the \(3D\times 2D\) area surrounding the skirted foundation, the mesh size was adjusted to \(0.5\times 0.5\) m to obtain more precise bearing capacity results. A total of 5300 four-node bilinear reduced integration  
Figure 1: A typical series of the soil random field: **(a)** Case Ani-1; **(b)** Case Ani-3; **(c)** Case Ani-6; **(d)** Case Ani-11. Figure 2: Details of skirted foundation and positive loading and displacement. quadrilateral (CPE4R) elements were used to simulate the seabed soil domain in the model, as shown in Figure 3.  
The soil was defined as an elastoplastic material that obeys the Tresca yield criterion, in which the maximum shear stress was determined by the undrained shear strength [27]. The soil undrained Young's modulus (\(E_{u}\)) was assumed to be linearly correlated to \(S_{u}\) with the relationship expressed in Equation (3):  
\[E_{u}=K_{u}\times S_{u} \tag{3}\]  
where \(K_{u}\) is the coefficient. In fact, the ultimate bearing capacity of foundation was confirmed to be barely influenced when \(K_{u}\) was equal to 500 [50]. In addition, the Poisson's ratio was fixed at 0.49 to simulate the undrained condition. As only the ultimate bearing capacities were of interest, the skirted foundation was simulated as a rigid body without any deformation. The skirted foundation and surrounding soil were coupled by the \"tie\" constraint to avoid the generation of extra interface element. The bearing capacity factors were introduced in this paper to evaluate the uniaxial bearing capacity normalized by the foundation diameter and the mean undrained shear strength at mudline, as expressed in Equation (4):  
\[ \begin{array}{l}N_{cH}=H_{0}/D\mu_{S_{u,0}}\\ N_{cV}=V_{0}/D\mu_{S_{u,0}}\\ N_{cM}=M_{0}/D^{2}\mu_{S_{u,0}}\\ \end{array}