## 2 Constitutive Model and Calibration  
### Constitutive Model  
The function of yield surface, considering the Von-Mises criterion, is defined by Equation (1) [30]:  
\[F=f\big{(}\sigma_{ij},\alpha_{ij}\big{)}=\sqrt{3J_{2}\big{(}\sigma_{ij}-\alpha _{ij}\big{)}}-\sigma_{0} \tag{1}\]  
where \(\sigma_{o}\) is the initial size of the yield surface under zero plastic strain, i.e., at the elastic stress limit, \(\alpha_{ij}\) is the back stress that refers to the center position of the yield surface, and \(J_{2}\big{(}\sigma_{ij}\big{)}\) is the second deviatoric stress invariant, given by Equation (2):  
\[J_{2}\big{(}\sigma_{ij}-\alpha_{ij}\big{)}=\frac{1}{2}(S_{mn}-t_{mn}).(S_{mn}- t_{mn}) \tag{2}\]  
where \(\sigma_{ij}\) and \(\alpha_{ij}\) are the stress and back stress tensor, respectively, and \(S_{mn}\) and \(t_{mn}\) are the stress and back stress deviatoric tensor, respectively. At an infinitely large plastic strain, when \(\sigma_{ij}\) approaches the effective maximum yield stress (\(\sigma_{y}\)) (Equation (3)), the back stress tensor (\(\alpha_{ij}\)) is given by Equation (4), as follows:  
\[\sigma_{y}=\sqrt{J_{2}\big{(}\sigma_{ij}\big{)}} \tag{3}\]  
\[\alpha_{ij}=\alpha_{s}=\sqrt{\frac{1}{2}(t_{mn}.t_{mn})}. \tag{4}\]  
The important aspect of the kinematic hardening model is the back stress evolution law, which describes the translation of the yield surface in the stress space. The kinematic hardening model consists of two parts--a purely kinematic term and a relaxation term--and can be written as (Equation (5)):  
\[d\alpha_{ij}=C\frac{1}{\sigma_{o}}\big{(}\sigma_{ij}-\alpha_{ij}\big{)}de^{pl }-\gamma_{o}\alpha_{ij}de_{pl} \tag{5}\]  
where \(C\) and \(\gamma_{o}\) are material parameters, \(de^{pl}\) is the plastic strain rate tensor, and \(de_{pl}\) is the equivalent plastic strain rate. In relaxation terms, the back stress is related to the equivalent plastic strain and the linear kinematic model is transformed to the nonlinear kinematic model. The plastic flow rule associated with the yield function is as follows (Equation (6)) [30]:  
\[de^{pl}=d\lambda\frac{\partial F}{\partial\sigma_{ij}} \tag{6}\]where the constant \(d\lambda\) is defined as \(d\lambda=d\varepsilon_{pl}\). The detailed expression can be found in Appendix A.  
Parameter \(\sigma_{0}\) represents the initial size of the yield surface at zero plastic strain and controls the onset of nonlinear soil behavior; it is defined as the fraction (\(\lambda\)) of the (\(\sigma_{y}\)) effective maximum yield stress at large plastic strain. Decreasing this ratio improves the model nonlinearity, but if it is too small, the model will be non-economic. In contrast, an excessive increase in this quantity has a negative effect on the model that leads to unrealistic results. In the 1D condition, it is given by Equation (7):  
\[\sigma_{0}=\lambda\sigma_{y} \tag{7}\]  
where \(\lambda\) is the material constant. The evolution Equation of kinematic hardening can be integrated to give Equation (8):  
\[\alpha=\frac{C}{\gamma_{o}}\Big{[}1-exp\Big{(}-\gamma_{o}\varepsilon_{pl} \Big{)}\Big{]} \tag{8}\]  
where \(\alpha\) is the back stress, \(\varepsilon_{pl}\) is the equivalent plastic strain, \(C\) is the initial kinematic hardening modulus (\(C=E=2(1+\
u)G_{0}\)) that defines the maximum translation of the yield surface, and \(\gamma_{o}\) defines the rate decrease of the kinematic hardening modulus. Parameter \(C\) can be derived by differentiating Equation (8) to give Equation (9):  
\[\frac{\partial\alpha}{\partial\varepsilon_{pl}}|_{\varepsilon_{pl}\to 0}C=E. \tag{9}\]  
The occurrence of displacement of the yield surface and the creation of back stress (\(\alpha_{ij}\)) in cyclic loading can be understood by considering the kinematic hardening model. At infinite large plastic strain (when \(\alpha_{ij}\) approaches \(\sigma_{y}\)), the back stress (\(\alpha_{ij}\)) is equal to \(\alpha_{s}=\frac{C}{\gamma_{o}}\) and (\(\sigma_{ij}-\alpha_{ij}\)) tends to \(\sigma_{0}\) (see Figure 1), which means that Equation (5) tends to zero.