## 2 Non-Equilibrium Thermodynamic-Based Model Framework
### Expression of Stress and Strain  
The elastic potential energy density function \(w_{e}\) is generally written as a function of elastic strain, such as the linear elastic model (generalized Hooke's law). In addition, according to the influence of average effective stress on elastic bulk modulus \(K_{e}\) and elastic shear modulus \(G_{e}\), some authors [29; 36] constructed the elastic potential energy density functions, which have similar forms to the linear elastic model. By analogy with Hertz contact, [PERSON] and [PERSON] [29] obtained the elastic potential energy density function of granular materials. \[w_{e}=B(\varepsilon_{v}^{e})^{0.5}\left[\frac{2}{5}(\varepsilon_{v}^{e})^{2}+ \frac{1}{\zeta}(\varepsilon_{s}^{e})^{2}\right] \tag{11}\]  
where \(\varepsilon_{v}^{e}\) is the elastic volume strain, \(\varepsilon_{s}^{e}\) is the elastic shear strain, \(\zeta\) is the material parameters, \(\zeta\) = 5/3 for dry sand, and \(B\) describes the hardness of the material, the same dimension as the stress. For unsaturated geotechnical materials, the cementation and the matrix suction has a great influence on the contact between particles, and cohesion between particles affects the elastic shear modulus. Following the method of describing the thermal expansion effect of solid [30], adding the contribution of suction \(\chi se_{v}^{e}\) into Equation (11), and considering the effect of cohesion, Equation (11) can be modified as follows:  
\[w_{e}=B\left[\frac{2}{5}(\varepsilon_{v}^{e})^{2.5}+\xi(\varepsilon_{v}^{e}+c )^{0.5}(\varepsilon_{s}^{e})^{2}\right]+\chi s\varepsilon_{v}^{e} \tag{12}\]  
where \(\xi=\)1/\(\zeta\); \(c\) is the material parameter to reflect the influence of cohesion on the elastic shear modulus and \(\chi=s_{r}\), which is similar to [PERSON]'s effective stress of unsaturated soil. \(\pi_{ij}\) is regarded as the mean effective stress of unsaturated geotechnical materials and is related to the elastic potential energy density \(w_{e}\). According to Equation (1), there is  
\[\pi_{ij}=\frac{\partial w_{e}}{\partial\varepsilon_{ij}^{e}}=B\left[( \varepsilon_{v}^{e})^{1.5}\delta_{ij}+0.5\xi(\varepsilon_{v}^{e}+c)^{-0.5}( \varepsilon_{s}^{e})^{2}\delta_{ij}+2\xi(\varepsilon_{v}^{e})^{0.5}\varepsilon _{ij}^{e}\right]+\text{ss}_{r}\delta_{ij} \tag{13}\]  
Similar to the linear elastic model, \(\pi_{ij}\) can be expressed as:  
\[\pi_{ij}=(K_{e}\varepsilon_{v}^{e}+\text{ss}_{r})\delta_{ij}+2G_{e}\varepsilon_ {ij}^{e} \tag{14}\]  
where the \(K_{e}\) and \(G_{e}\) can be expressed as:  
\[K_{e}=B\left[(\varepsilon_{v}^{e})^{0.5}+0.5\xi(\varepsilon_{v}^{e}+c)^{-0.5} \frac{(\varepsilon_{s}^{e})^{2}}{\varepsilon_{v}^{e}}\right] \tag{15}\]  
\[G_{e}=B\xi(\varepsilon_{v}^{e})^{0.5} \tag{16}\]Additionally, the mean effective stress \(p^{\prime}\) and shear stress \(q\) can be obtained:  
\[p^{\prime}=\frac{\pi_{kk}}{3} =B\Big{[}\big{(}\epsilon_{v}^{e})^{1.5}+0.5\xi(\epsilon_{v}^{e}+c)^ {-0.5}(\epsilon_{s}^{e})^{2}\Big{]}+\text{ss}_{r} \tag{17}\] \[q =\sqrt{\frac{3}{2}s_{ij}s_{ji}}=\sqrt{6}B\xi(\epsilon_{v}^{e}+c)^ {0.5}\epsilon_{s}^{e} \tag{18}\]  
where \(s_{ij}=\pi_{ij}-\frac{\pi_{kk}}{3}\delta_{ij}\) is the deviator stress. It can be seen that \(B\) is related to the density, and the greater the density, the greater the \(B\). Thus, \(B=B_{0}\exp\left(-\frac{c}{\lambda(s)}\right)\) can be obtained according to the compression curve (\(e\)-\(\text{ln}p^{\prime}\)) of geotechnical materials [26], where \(e\) is the void ratio and \(\lambda(s)\) is the slope of compression curve. According to the migration coefficient matrix, the following relationship can be obtained:  
\[\dot{\epsilon}_{ij}^{p}=\eta_{ijkl}\pi_{kl} \tag{19}\]  
When the granular temperature \(T_{g}=0\), the plastic strain will not occur.