## 2 Constitutive Model for Undrained Clay
### The Non-Iterative Von Mises Model of [PERSON] and [PERSON] [39]  
Stress and strain tensors are decomposed into volumetric and deviatoric parts, i.e.,  
\[\mathbf{\sigma}=\mathbf{s}+p\mathbf{I} \tag{1}\]  
\[\mathbf{\varepsilon}=\mathbf{\varepsilon}_{d}+\frac{1}{3}\mathbf{\varepsilon}_{v}\mathbf{I} \tag{2}\]  
where \(p\) and \(\mathbf{\varepsilon}_{v}\) are the mean stress and volumetric strain, \(\mathbf{s}\) and \(\mathbf{\varepsilon}_{d}\) are the second-order tensors of deviatoric stress and shear strain, and \(\mathbf{I}\) is the second-order unit tensor (tensor quantities are shown by boldface font, while scalars are shown by regular font). The volumetric and shear deformations are given by  
\[p=K\mathbf{\varepsilon}_{v} \tag{3}\]  
\[\mathbf{s}=2G\Big{[}\mathbf{\varepsilon}_{d}-\mathbf{\varepsilon}_{d}^{p}\Big{]} \tag{4}\]  
where \(K\) and \(G\) are the bulk and shear moduli, respectively, and \(\mathbf{\varepsilon}_{d}^{p}\) is the second-order tensor of plastic shear strains. In order to have kinematic hardening, a second-order tensor of back-stress (\(\mathbf{\alpha}\)) needs to be defined as the center of the von Mises yield surface in the deviatoric stress space. A relative stress tensor is defined as \(\mathbf{\Sigma}=\mathbf{s}-\mathbf{\alpha}\). The von Mises yield criterion can be defined by  
\[f=\|\mathbf{\Sigma}\|-\sqrt{\frac{2}{3}}\mathbf{\sigma}_{y} \tag{5}\]  
where \(\|\mathbf{\Sigma}\|\) is the norm of \(\mathbf{\Sigma}\), \(\mathbf{\sigma}_{y}=\mathbf{\sigma}_{y0}+H_{iso}\mathbf{\varepsilon}_{d}^{p}\), \(\mathbf{\sigma}_{y0}\) is initial yield stress, \(H_{iso}\) is the isotropic hardening coefficient, and \(\mathbf{\varepsilon}_{d}^{p}\) is the equivalent accumulated plastic strain (EAPS) defined by  
\[\mathbf{\varepsilon}_{d}^{p}\equiv\int_{0}^{t}\sqrt{\frac{2}{3}}\Big{\|}\mathbf{ \varepsilon}_{d}^{p}\Big{\|}dt \tag{6}\]  
where \(\mathbf{\dot{\varepsilon}}_{d}^{p}\) is the rate of the plastic shear strains tensor, which, with the assumption of associative plasticity, can be obtained by \(\mathbf{\dot{\varepsilon}}_{d}^{p}=\dot{\gamma}\mathbf{n}\), with \(\dot{\gamma}\) being the plastic rate that can be related to EAPS by \(\mathbf{\dot{\varepsilon}}_{d}^{p}=\sqrt{\frac{2}{3}}\dot{\gamma}\) and \(\mathbf{n}\) being the second-order tensor defined as \(\mathbf{n}\equiv\frac{\mathbf{\pi}}{\|\mathbf{\Sigma}\|}\). Kinematic hardening is defined by the evolution of the back-stress, i.e.,  
\[\dot{\mathbf{\alpha}}=\frac{2}{3}H_{kin}\mathbf{\dot{\varepsilon}}_{d}^{p}-H_{nl}\mathbf{ \dot{\varepsilon}}_{d}^{p}\mathbf{\varepsilon} \tag{7}\]  
where \(\dot{\mathbf{\alpha}}\) is the rate of the back-stress tensor, and \(H_{kin}\) and \(H_{nl}\) are the kinematic and non-linear hardening coefficients, respectively. For a given strain, first a trial stress is calculated by Equations (3) and (4), assuming the \(\epsilon_{d}^{p}\) retains its value from the previous time increment. Then, the yield function is calculated with Equation (5), maintaining the value of EAPS and back-stress from the previous increment. If the yield function is negative, the assumption of elasticity is confirmed. Otherwise, return-mapping is carried out as follows. The discrete form of a plastic multiplier (\(\lambda=\hat{\gamma}dt\)) is obtained by finding the smallest positive real root of the fourth-degree algebraic equation, i.e.,  
\[C_{1}\lambda^{4}+C_{2}\lambda^{3}+C_{3}\lambda^{2}+C_{4}\lambda+C_{5}=0 \tag{8}\]  
The coefficients \(C_{1}\)-\(C_{5}\) are functions of the shear modulus \(G\), hardening moduli \(H_{iso}\), \(H_{kin}\) and \(H_{nl}\), yield stress from the previous increment \(\sigma_{ym}\), and trial tensors of deviatoric stress \(\mathbf{s}^{TR}\) and back-stress \(\mathbf{a}^{TR}\). The functions for the coefficients are given in [39] and are not repeated here for brevity.