## 2 Numerical Methodology
### Constitutive Models  
For drained shear of dilating sands, it is common practice to compute the peak friction angle \(\phi_{P}\) as the sum of the critical state friction angle \(\phi_{C}\) and the dilatancy term \(\psi^{\prime}\) which depends on void ratio \(e\) and effective mean pressure \(p^{\prime}\).  
\[\phi^{\prime}=\phi_{c}+\psi^{\prime}\big{(}p^{\prime},e\big{)} \tag{2}\]  
A widely used expression in the form of Equation (2) is that of [PERSON] [36], which can be put in the form:  
\[\phi_{P}=\phi_{C}+\triangle\phi I_{R}=\phi_{C}+\triangle\phi\Big{(}I_{D}\Big{(} Q-\ln\Big{(}p^{\prime}/p^{\prime}_{ref}\Big{)}\Big{)}-1\Big{)},\ 0\leq I_{R}\leq 4 \tag{3}\]  
Using the data of 17 soils, [PERSON] [36] proposed that the fit parameter \(Q=10\), \(\triangle\phi=5^{\circ}\) for plane strain and \(\triangle\phi=3^{\circ}\) for triaxial strain, with reference pressure \(p^{\prime}_{ref}\) taken as 1 kPa. And the peak dilation angle \(\psi_{P}\) is also a function of the relative dilatancy index \(I_{R}\):  
\[\psi_{P}=0.8\bigtriangleup\phi I_{R}=0.8\bigtriangleup\phi\Big{(}I_{D}\Big{(} Q-\ln\Big{(}p^{\prime}/p^{\prime}_{ref}\Big{)}\Big{)}-1\Big{)} \tag{4}\]  
in which \(I_{D}\) is the relative density of sand. Negative value of \(I_{R}\) indicates that large contractile strains will occur before \(\phi_{C}\) can be mobilized, which typically happens in loose sand.  
For positive value of \(I_{R}\), the effective friction angle \(\phi^{\prime}\) is larger than the critical friction angle \(\phi_{C}\). However, if progressive failure occurs, all \(\phi^{\prime}\) values in excess of \(\phi_{C}\) are transitory with larger strains, which is defined as strain softening behavior of the sand. As depicted in Figure 3, a simple model proposed by [PERSON] et al. [37] is adopted to describe the softening behavior.  
In the softening model shown in Figure 3, the effective friction angle \(\phi^{\prime}\) and dilation angle \(\Psi^{\prime}\) changes with the accumulated deviatoric plastic strain \(k_{shear}\). \(\phi^{\prime}\) decreases from \(\phi_{p}\) to \(\phi_{c}\) as \(k_{shear}\) increases from \(k_{shear}^{p}\) to \(k_{shear}^{c}\), whereas \(\Psi^{\prime}\) decreases from \(\Psi_{p}\) to 0. The parameter \(k_{shear}\) is defined as follows:  
\[k_{shear}=\int\dot{k}_{shear}dt \tag{5}\]  
in which  
\[\dot{k}_{shear}=\sqrt{0.5e_{ij}^{p}e_{ij}^{p}} \tag{6}\]  
\(t\) is time, and \(\dot{k}_{shear}\) is the deviatoric plastic strain rate tensor, of which the expression is  
\[e_{ij}^{p}=\varepsilon_{ij}^{p}-\frac{1}{3}\varepsilon_{kk}^{p}\dot{\delta}_ {ij} \tag{7}\]  
where \(\delta_{ij}\) is the Kronecker tensor, and \(\varepsilon_{kk}^{p}\) is the plastic deformation rate tensor.  
Spudcan penetration involves large plastic strain in clay, which means that the strain softening effect of clay should be considered in the present model. The model proposed by [PERSON] and [PERSON] [28] is incorporated into the present model as follows:  
\[S_{u}=\bigg{[}\delta_{rem}+(1-\delta_{rem})\exp\left(-\frac{3\xi}{\xi_{95}} \right)\bigg{]}S_{u0} \tag{8}\]  
in which \(S_{u}\) is the undrained shear strength of clay used for simulation, while \(S_{u0}\) is the initial undrained shear strength of the clay; \(\delta_{rem}\) is the residual degradation factor which is the inverse of soil sensitivity St; \(\xi\) is the accumulated plastic shear strain; and \(\xi_{95}\) is the reference accumulated plastic shear strain which is taken as 10 [38; 39].