## 2 Model Description
Then the deviatoric and volumetric strain increment can be composed as follows:  
\[\dot{\epsilon}_{\mathrm{q}}=\dot{\epsilon}_{\mathrm{q}}^{\mathrm{e}}+\dot{ \epsilon}_{\mathrm{q}}^{\mathrm{P}}=\frac{\mathrm{q}\cdot}{3\mathrm{G}}+\frac{ \mathrm{p}^{\prime}\cdot\dot{\eta}}{\mathrm{K}_{\mathrm{P}}} \tag{9}\]  
\[\dot{\epsilon}_{\mathrm{v}}=\dot{\epsilon}_{\mathrm{v}}^{\mathrm{e}}+\dot{ \epsilon}_{\mathrm{v}}^{\mathrm{P}}=\frac{\dot{\mathrm{p}^{\prime}}}{\mathrm{K }}+\mathrm{d}\frac{\mathrm{p}^{\prime}\cdot\dot{\eta}}{\mathrm{K}_{\mathrm{P}}} \tag{10}\]  
Here, it is of great importance to recall that the flow rule used in this study is non-associated, as the flux of plastic strain is directly defined by the dilatancy equation. Meanwhile, the flux of plastic strain is controlled by model constant parameter \(\mathrm{n}\). State parameter is implanted into plastic hardening modulus \(\mathrm{K}_{\mathrm{P}}\)[9; 12]:  
\[\mathrm{K}_{\mathrm{P}}=\mathrm{CG}\Big{(}\frac{\mathrm{M}}{\mathrm{n}}-\psi \Big{)} \tag{11}\]  
where \(\mathrm{G}\) is the elastic shear modulus at the current state. Note that the initial state parameter values \(\psi_{\mathrm{ini}}=\mathrm{e}_{0}/\mathrm{e}_{c}\), under different density index and confining pressures are possible to them to be the same (see Figure 1). Thus, a dimensionless model parameter \(C(\mathrm{e}_{0})\) standing for the effect of initial void ratio is needed in Equation (11) in order to compensate this error, and it reads:  
\[C=\alpha\mathrm{e}_{0}^{\beta} \tag{12}\]  
where \(\alpha\), \(\beta\) are model dimensionless parameters and \(\mathrm{e}_{0}\) is the initial void ratio. It is worth noting that the plastic hardening modulus \(K_{\mathrm{P}}\) has the following special values:  
* \(K_{\mathrm{P}}(\eta=0)=\lim\limits_{\eta=0}\bigl{(}\mathrm{C}G\bigl{(}\frac{M}{ \eta}-\psi\bigr{)}\bigr{)}=\infty\), at the beginning of the test;
* \(K_{\mathrm{P}}(\eta=M,\psi=1)=\lim\limits_{\eta=M;\psi=1}\bigl{(}\mathrm{C}G \bigl{(}\frac{M}{\eta}-\psi\bigr{)}\bigr{)}\), in achievement of the critical state. From Equations (11) and (12), it can be seen that the hardening or softening behaviour is determined by the sign of \(\bigl{(}\frac{M}{\eta}-\psi\bigr{)}\). Sand is initially in a dense state when \(\psi<1\), which results in \(\frac{M}{\eta}-\psi>0\) and sand response is hardening. Sand is initially in a loose state when \(\psi<1\). Its response is hardening with \(\bigl{(}\frac{M}{\eta}-\psi\bigr{)}>0\) and softening with \(\bigl{(}\frac{M}{\eta}-\psi\bigr{)}<0\). Thus, the introduction of state parameter in the plastic modulus succeeds in distinguishing the hardening and softening response of sand. Finally, the relationship between the stress and strain increments can be established in a matrix form as:  
\[\left( \begin{array}{c}\dot{\epsilon}_{\mathrm{q}}\\ \dot{\epsilon}_{\mathrm{v}}\end{array} \right)=\left[S\right]\times\left( \begin{array}{c}\dot{\mathrm{q}}\\ \dot{\mathrm{P}}\end{array} \right)=\left[ \begin{array}{cc}\frac{1}{3C}+\frac {1}{K_{\mathrm{P}}}&-\frac{\eta}{K_{\mathrm{P}}}\\ \frac{d}{K_{\mathrm{P}}}&\frac{1}{K}-\frac{d\eta}{K_{\mathrm{P}}}\end{array}  \right]\times\left( \begin{array}{c}\dot{\mathrm{q}}\\ \dot{\mathrm{P}}\end{array} \right) \tag{13}\]  
\[\left( \begin{array}{c}\dot{\mathrm{q}}\\ \dot{\mathrm{P}}\end{array} \right)=\left[S\right]^{-1}\times\left( \begin{array} []{c}\dot{\epsilon}_{\mathrm{q}}\\ \dot{\epsilon}_{\mathrm{v}}\end{array} \right)=\left[ \begin{array}{cc}3G+9 AG^ {2}&-3 AKG\eta\\ 3 AKG\mathrm{d}&K-AK^{2}\eta\end{array} \right]\times\left( \begin{array}{c} \dot{\epsilon}_{\mathrm{q}}\\ \dot{\epsilon}_{\mathrm{v}}\end{array} \right) \tag{14}\]  
where A is equal to \(1/\bigl{(}K\eta\mathrm{d}-3G-K_{\mathrm{P}}\bigr{)}\).