1 Introduction

It has widely acknowledged that stress–strain behaviour of geomaterials (such as rockfills, rocks, and soils) is affected by loading time or rate—in other words geomaterials are time-dependent materials in general [13, 14, 17, 22, 29, 43, 47, 50, 58, 60, 73].

Time-dependent mechanical behaviour of soils in the saturated state has been widely reported, particularly soft clayey soils. For example, many experimental studies [4, 9, 19, 48, 55, 56, 71] have been performed to investigate the creep behaviour of saturated soils under controlled stresses in laboratory. In addition to creep, time-dependent behaviour of soils can also be reflected through investigating their behaviour with different strain rates [18, 25, 32, 44, 52, 66, 72]. For instance, it is found that there is an increase in the apparent preconsolidation pressure of saturated soils when strain rates increase [18, 25, 33]. Based on experimental results and mechanism analyses, numerous constitutive models have been developed to mimic and explain the observed time-dependent mechanical behaviour of saturated soils since 1960s [1, 5, 10, 16, 20,21,22,23, 26, 30, 39, 45, 47, 54, 60, 65, 67]. Herein, the most fundamental work was pioneered by Perzyna [39] who proposed perhaps the first elasto-viscoplastic (EVP) model to explain the time-dependent behaviour for saturated soils. Various EVP constitutive models have been proposed on the basis of overstress theory from Perzyna’s study [39], such as Adachi and Oka [1], Borja and Kavazanjian [10], Fodil et al. [16], Madaschi and Gajo [30], Yin et al. [67] and Sun and Sumelka [47]. Following the experimental observation that loading lines with different loading rates are almost parallel to each other, Bjerrum [11] has proposed the concept of equivalent (or reference) time, and a large number of EVP constitutive models for saturated soils have been formulated by using this concept [10, 20,21,22,23, 54, 60, 65, 67].

In addition to saturated soils, in fact, the mechanical behaviour of unsaturated soils is also found time-dependent [12, 13, 17, 27, 41, 50, 57, 68]. Many researchers [2, 15, 24, 27, 42, 64] investigated creep characteristics of unsaturated soils through experimental studies. They have quantified the dependency of the viscous parameter on suction, and proposed creep models for unsaturated soils are at different suctions based on the experimental measurements. Similar to saturated soils, mechanical behaviour of unsaturated soils also shows strain rate dependency (or time-dependency in general) according to the experimental studies in the literature [6, 7, 17, 36, 41, 50, 57]. For example, one-dimensional oedometer tests done by Qin et al. [41] have showed that, under the same suction (s), loading curves with different strain rates run log-linear and parallel and yield stress of unsaturated soils increases with strain rates, as observed by other researchers [17, 38]. This observation implies that the isotach concept for saturated soils can also be applied to unsaturated soils. Wu et al. [57] studied the hydro-mechanical behaviour of a reconstituted unsaturated soil under different suctions and strain rates through various rate-controlled unsaturated/undrained triaxial tests. They concluded that, like saturated soils, the critical state lines for unsaturated soils with constant strain rates are parallel with each other in the e − lnp space. On the basis of the experimental observations, Gennaro and Pereira [17] have made attempts to study time-dependent stress–strain behaviour of unsaturated geomaterials with EVP constitutive frameworks. The Barcelona Basic Model (BBM) from Alonso et al. [3] has been extended to an EVP model for unsaturated chalks with the hardening law including rate effects. It should be noted that the EVP model developed by Gennaro and Pereira [17] is restricted to normally consolidated unsaturated soils, and loading time histories are overlooked in the model. However, many experimental and constitutive studies have demonstrated that stress histories significantly affect behaviour of unsaturated soils [31, 34, 35, 49, 51, 53, 61, 70]. Besides, the formulation to account for the dependency of the viscous parameter on suction in their EVP framework still have room to be improved to reflect the observations from more recent in experimental data [2, 15, 17, 41, 64].

The literature review suggests that, numerous EVP constitutive models have been reported with regarding to saturated soils, but there are very limited studies focusing on the EVP constitutive frameworks for unsaturated soils. The EVP constitutive model for unsaturated soils proposed by Gennaro and Pereira [17] is a pioneering work but still suffers from some drawbacks as mentioned above. In practice, due to complex environments, unsaturated soils inevitably undergo different loading and creep histories and are subjected to different natural or man-induced saturation-desaturation processes (such as quick wetting due to a heavy rainstorm and slow wetting due to long-last small rain). Consequently, how loading and creep histories, and drying/wetting rates affect EVP behaviour of unsaturated soils needs to be carefully assessed. To better address these issues and overcome the limits of existing models, an advanced constitutive framework for unsaturated soils, the UTUH (Unsaturated Time-dependent Unified Hardening) model, is proposed in this paper to consider the joint effect of time, suction and overconsolidation based on sub-loading surface plasticity, EVP framework and unsaturated soil mechanics. This paper is organized as follows. First, a reference line, namely the instantaneously normal compression line (INCLs) for unsaturated soils, is introduced to define creep time and overconsolidation states of unsaturated soils. Subsequently, an isotropic EVP constitutive framework for unsaturated soils is proposed by combining viscous deformation with mechanical and hydraulic deformation through overconsolidation parameter. In this paper, net stress, suction and time are adopted as fundamental constitutive variables because the stress-suction approach not only contributes to capturing typical behaviour of unsaturated soils, but is consistent with common laboratory tests [70], and time-dependent LC yield surface is derived to characterize the relationship between yield stress, suction and time. We extend the time-dependent LC yield surface to triaxial stress state and involve sub-loading surface to assess the effect of overconsolidation. The hardening of yield surface and sub-loading surface is controlled by viscoplastic volumetric strain and unified hardening parameter [59, 60, 62]. The UTUH model, which only requires 13 parameters, was calibrated based on the experimental data from the literature. Model performance discussions and systematic experimental validations are presented before the concluding remarks.

2 Isotropic EVP constitutive model for unsaturated soils

2.1 Time-dependent deformation

Qin et al. [41] have conducted a series of oedometer tests for GMZ01 bentonite under controlled suctions and CRSs (constant rates of strain). Three loading rates, 1 × 10–5, 1 × 10–6 and 1 × 10–7/s, were adopted, and the loading strategies were divided into two categories: the single-stage CRS and the stepwise CRS. The test results are shown in Fig. 1.

Fig. 1
figure 1

CRS test results under controlled suctions (data after Qin et al. [41])

Some key observations can be summarized from Fig. 1 as follows: (1) Under the same suction (s), loading lines with different loading rates run nearly parallel to each other. In other words, the isotach concept for saturated soils, e.g. a series of time lines in Yao et al. [60], can also be applied to unsaturated soils. This observation agrees with study results reported by other researchers [17, 38]. (2) Under constant rates, the slope for loading lines (i.e. Cc3, Cc2 and Cc1) decreases with the increase of suction, which is identical to the conclusion drawn by Alonso et al. [3] when rate is ignored. (3) The intervals between adjacent loading lines (constant strain rates) become narrow with the increase of suction—implying time-dependency becomes weak with increase of suction and coefficient of the secondary consolidation Cαe decreases when drying, as shown in Fig. 2. The decrease in viscous property was also observed and addressed by other researchers [13, 17, 38, 40, 68]. Qin et al. [41] suggested that Cαe can be expressed as Eq. (1):

$$ C_{{\alpha {\text{e}}}} \left( s \right) = C_{c} \left( s \right)\varphi \left( s \right) $$
(1)

where \(\varphi \left( s \right) = \varphi \left( 0 \right) - b\ln \left( {\frac{{s + p_{{{\text{at}}}} }}{{p_{{{\text{at}}}} }}} \right)\), φ(0) is ratio of coefficient of the secondary consolidation to compression coefficient Cc when suction is 0 (saturated soils), b is the parameter defining the effect of suction on Cαe, and pat denotes atmospheric pressures. For isotropic states, Eq. (1) can also be expressed as:

$$ \beta \left( s \right) = \lambda \left( s \right)\varphi \left( s \right) $$
(2)

where λ(s) and β(s) are elastoplastic compression index and the viscous parameter for unsaturated soils with suction s, respectively. Here we assume φ(0) and b can still be used for isotropic states.

Fig. 2
figure 2

The effect of suction on Cαe (data after Qin et al., [41])

According to the phenomena observed in Fig. 1, a conceptual framework (see Fig. 3) is proposed here for unsaturated soils following authors’ previous study on time-dependency of saturated soil [60]. Any laboratory test of unsaturated soils (including saturated soil as a special case) involves time. Therefore, not only mechanical deformation (including suction-induced compression like drying and wetting), but time-dependent deformation, as long as the mechanical loading takes time to complete. As shown in Fig. 3, void ratio theoretically follows the path \(\overline{\mathrm{AB} }\) during instant isotropic compression under suction s for unsaturated soils with zero loading time, but in fact the time effect (e.g. deformation \(\overline{\mathrm{BC} }\) in Fig. 3) leads to the path \(\overline{\mathrm{AC} }\) in reality. More generally, NCLs (i.e. NCL with constant suction s) obtained in the laboratory test inevitably involves time-dependent deformation. Therefore, when the loading rate is infinite or loading time is infinitesimal, NCLs becomes INCLs which is a theoretical line under instant compression condition, which represents a benchmark with zero time effect. INCLs can be expressed as

$$e=e\left(s\right)-\lambda \left(s\right)\mathrm{ln}\frac{p}{{p}^{\mathrm{c}}}$$
(3)

where e is void ratio, p is net mean stress, pc is reference stress, and e(s) is the void ratio at reference stress on the INCLs of unsaturated soils where \(e\left(s\right)=e\left(0\right)-{\kappa }_{\mathrm{s}}\mathrm{ln}\frac{s+{p}_{\mathrm{at}}}{{p}_{\mathrm{at}}}\) [3] where κs is an elastic stiffness parameter with respect to suction change.

Fig. 3
figure 3

The INCLs in a conceptual framework

To obtain time-dependent deformation, Eq. (4) is proposed here after previous study for saturated soils done by Yao et al. [60].

$$\Delta e=-\beta \left(s\right)\mathrm{ln}\left(\frac{{t}_{a}}{{t}_{0}}+1\right)$$
(4)

where t0 is reference time to eliminate the effect of unit of time, which is taken as 1 min in this paper [60] and ta is creep time involving creep history. Since for unsaturated soil, viscous parameter β is suction-dependent here. Under constant suction, two loading paths, the “pure” creep path (i.e. the path \(\overline{\mathrm{BD} }\) in Fig. 3) and the instant mechanical loading–unloading path (i.e. the path \(\overline{\mathrm{BED} }\) in Fig. 3) are chosen to calculate plastic deformation, respectively. Hence, there is

$$ \Delta e_{{\left| {{\text{BD}}} \right.}} = - \beta \left( s \right)\ln \left( {\frac{{t_{a} }}{{t_{0} }} + 1} \right) = \Delta e_{{\left| {{\text{BED}}} \right.}} = - \left[ {\lambda \left( s \right) - \kappa } \right]\ln \frac{{p_{E} }}{{p_{{{\text{st}}}} }} = - \left[ {\lambda \left( s \right) - \kappa } \right]\ln {\text{OCR}} = \left[ {\lambda \left( s \right) - \kappa } \right]\ln R $$
(5)

where pE and pst are net mean stress at the points E and D, respectively. OCR is the overconsolidation ratio in terms of mean stress, and κ is elastic compression index, which is assumed to be suction and time-independent here. R is the overconsolidation parameter and is defined by reciprocal of the overconsolidation ratio if deviator stress is equal to zero. According to Eq. (5), ta is expressed as:

$${t}_{a}={t}_{0}\left({R}^{-\mathrm{\alpha }}-1\right)$$
(6)

where \(\alpha =\frac{\lambda \left(s\right)-\kappa }{\beta \left(s\right)}\). According to Eqs. (4) and (6), time-dependent volumetric deformation can be written as:

$$ {\text{d}}\varepsilon_{v}^{t} = \frac{\beta \left( s \right)}{{1 + e_{0} }}\frac{{R^{\alpha } }}{{t_{0} }}{\text{d}}t $$
(7)

where e0 denotes initial void ratio, and t is creep time from instant states.

2.2 L-C and S-I yield curves

The BBM [3] proposed two yield curves, L-C curve and S-I curve, which defined the boundary between the elastic domain and the elastoplastic domain together. They are also adopted here, which can be expressed by Eqs. (8) and (9), respectively.

$$ \frac{{p_{s0} }}{{p^{{\text{c}}} }} = \left( {\frac{{p_{00} }}{{p^{{\text{c}}} }}} \right)^{\gamma } $$
(8)

where \(\gamma =\frac{\lambda \left(0\right)-\kappa }{\lambda \left(s\right)-\kappa }\),\(\lambda \left(s\right)=\lambda \left(0\right)\left[\left(1-r\right)\mathrm{exp}\left(-\xi s\right)+r\right]\), p00 and ps0 are consolidation pressure (i.e. yield mean stress) when soils are under saturated condition and unsaturated condition, respectively. λ(0) is elastoplastic compression index for saturated soils. r and ξ are two parameters defined in the BBM.

$$s={s}_{0}$$
(9)

where s0 is the yield suction of S-I yield curve.

2.3 Hardening laws

The hardening of L-C curve and S-I curve can be written as:

$$ \frac{{{\text{d}}p_{00} }}{{p_{00} }} = \frac{1}{{c_{p} }}{\text{d}}\varepsilon_{{\text{v}}}^{{{\text{tp}}}} $$
(10)
$$ \frac{{{\text{d}}s_{0} }}{{s_{0} + p_{{{\text{at}}}} }} = \frac{{1 + e_{0} }}{{\lambda_{{\text{s}}} - \kappa_{{\text{s}}} }}{\text{d}}\varepsilon_{{\text{v}}}^{{{\text{tp}}}} $$
(11)

where \({c}_{\mathrm{p}}=\frac{\lambda \left(0\right)-\kappa }{1+{e}_{0}}\), d\({\varepsilon }_{\mathrm{v}}^{\mathrm{tp}}\) is viscoplastic volumetric strain increment, and λs is an elastoplastic stiffness parameter with respect to suction change. From Eqs. (7), (10) and (11), it can be found that viscoplastic volumetric strain plays a key role in terms of hardening process. d\({\varepsilon }_{\mathrm{v}}^{\mathrm{tp}}\) can be produced due to mechanical loading, suction change, and creep. d\({\varepsilon }_{\mathrm{v}}^{\mathrm{tp}}=\mathrm{ d}{\varepsilon }_{\mathrm{v}}^{\mathrm{t}}\) if both stress and suction keep unchanged.

2.4 Elastic deformation

The increment for elastic volumetric strain, which is independent of time [65], can be written as Eq. (12):

$$\mathrm{d}{\varepsilon }_{\mathrm{v}}^{\mathrm{e}}=\frac{\kappa }{1+{e}_{0}}\frac{\mathrm{d}p}{p}+\frac{{\kappa }_{\mathrm{s}}}{1+{e}_{0}}\frac{\mathrm{d}s}{s+{p}_{\mathrm{at}}}$$
(12)

3 The UTUH model in triaxial stress state

To involve the overconsolidation effect which can be caused by stress/suction history and creep/time-dependency in triaxial stress state, a yield surface and a sub-loading surface are proposed here. The former represents the normally consolidated state, and the latter describes the overconsolidation state. The yield surface and the sub-loading surface dynamically evolve with the loading or unloading process.

3.1 Sub-loading surface

For saturated soil, Yao et al. [60] have developed a sub-loading surface to consider joint mechanical and rheological loading/unloading. Modification is made to combine sub-loading surface for saturated soil with Eq. (8) to include the suction effect following the BBM [3]. Regarding the hardening law for sub-loading surface, previous study has demonstrated that the unified hardening (UH) parameter H proposed by authors [59] can be employed to describe the hardening of sub-loading surface of soils subjected to time-dependency [60], temperature change [62], unsaturated state [70],and frozen conditions [46]. To consider joint effect of rheological load (i.e. time), hydraulic load (i.e. suction) and mechanical load (i.e. stress), UH parameter H is adopted to control the hardening of proposed sub-loading surface:

$$ \left\{ {\begin{array}{*{20}c} {f = {\rm ln}\left( {\frac{{q^{2} }}{{M^{2} p^{\prime\prime}}} + p} \right) + \frac{{\tilde{t}}}{\alpha } - {\rm ln}p_{{{\text{s}}0}} = 0} \\ {\frac{{p_{{{\text{s}}0}} }}{{p^{{\text{c}}} }} = \left( {\frac{{p_{00} }}{{p^{{\text{c}}} }}} \right)^{\gamma } } \\ {c_{p} {\text{d}}\left( {{\rm ln}p_{00} } \right) = {\text{d}}H} \\ \end{array} } \right. $$
(13)

where

$$ p^{\prime\prime} = p + p_{{\text{s}}} $$
$$ p_{{\text{s}}} = ks $$
$$ H = {\text{ }}\int {\frac{{M_{{\text{f}}}^{4} - \eta ^{4} }}{{M^{4} - \eta ^{4} }}{\text{d}}\varepsilon _{{\text{v}}}^{{{\text{tp}}}} } $$
$$ M_{{\text{f}}} \,{ = }\,6\left[ {\sqrt {\frac{\chi }{R}\left( {1{ + }\frac{\chi }{R}} \right)} - \frac{\chi }{R}} \right] $$
$$ \chi = \frac{{M^{2} }}{{12\left( {3 - M} \right)}} $$

where q is deviator stress, k is a parameter, and η is stress ratio \( \left( {\eta = \frac{q}{{p\prime \prime }}} \right) \). M is critical state stress ratio, and according to test results [7, 8, 69], it is assumed to be suction and time-independent here. Mf denotes potential failure stress ratio, and it is time-independent [60] and suction-independent [61]. \(\widetilde{t}\) is a time factor proposed by Yao et al. [60] initially for saturated soils to represent the contribution of time to the hardening process. An isotropic creep path can be chosen to quantify its expression. Hence, with Eqs. (7) and (13), there is:

$$\frac{1}{\alpha }\widetilde{t}=\int \frac{1+{e}_{0}}{\lambda \left(s\right)-\kappa }\mathrm{d}H$$
(14)
$$\widetilde{t}=\int \frac{{M}_{\mathrm{f}}^{4}}{{M}^{4}}\frac{{R}^{\mathrm{\alpha }}}{{t}_{0}}\mathrm{d}t$$
(15)

From Eq. (15), we can find the time factor is dimensionless variable related to time (t) and similarity ratio (R). Because the viscous deformation can be different at the same given time if soil’s state is different (such as different stress histories), time factor has a multiplier related to stress history (i.e. variable R which is the similarity ratio between sub-loading surface and yield surface to represent the relationship between current stress state and memorized stress state in the history). Specially, \(\widetilde{t}\) = \(\frac{t}{{t}_{0}}\) if soil is normally consolidated. If we assume that \(\frac{{q}^{2}}{{M}^{2}{p}^{{\prime}{\prime}}}+p\) is equivalent to pst (see Fig. 3 for isotropic stress state), Eq. (13) can be expressed as:

$${p}_{\mathrm{st}}={p}_{\mathrm{s}0}\mathrm{exp}\left(-\frac{\widetilde{t}}{\alpha }\right)={p}^{\mathrm{c}}{\left(\frac{{p}_{00}}{{p}^{\mathrm{c}}}\right)}^{\gamma }\mathrm{exp}\left(-\frac{\widetilde{t}}{\alpha }\right)$$
(16)

Equation (16) represents how pst evolves with time and suction—a time-dependent LC surface used as sub-loading surface in space of p, \(\widetilde{t}\) and s (see Fig. 4).

Fig. 4
figure 4

Sub-loading surface and yield surface in space of p, \(\widetilde{t}\) and s

3.2 Yield surface

Yield surface represents the normally consolidated state. The previous study [60] on saturated soil’s time-dependent behaviour suggested that the yield surface used in the Cam-clay model can be directly employed as the yield surface in the sub-loading elastoplastic model. Here, since unsaturated soil’s time-dependent behaviour is concerned, we employed the yield surface of the BBM model as the yield surface here, only updating plastic volumetric strain as viscoplastic volumetric strain to involve time-dependency:

$$ \left\{ {\begin{array}{*{20}c} {\bar{f} = {\rm ln}\left( {\frac{{\bar{q}^{2} }}{{M^{2} \bar{p}^{{\prime \prime }} }} + \bar{p}} \right) - {\rm ln}\bar{p}_{{{\text{s}}0}} = 0} \\ {\bar{p}_{{{\text{st}}}} = \bar{p}_{{{\text{s}}0}} = p^{{\text{c}}} \left( {\frac{{\bar{p}_{{00}} }}{{p^{{\text{c}}} }}} \right)^{\gamma } } \\ {c_{p} d\left( {{\rm ln}\bar{p}_{{00}} } \right) = d\varepsilon _{{\text{v}}}^{{{\text{tp}}}} } \\ \end{array} } \right. $$
(17)

where the over bar ‘-’ is used here to distinguish stress variables for the yield surface from that for the sub-loading surface. It should be noted that, since yield surface does not explicitly contain the time factor \(\widetilde{t}\) as in the sub-loading surface, the time-dependent behaviour has been involved in \(\mathrm{d}{\varepsilon }_{\mathrm{v}}^{\mathrm{tp}}\). Once both yield surface and sub-loading surface are defined, the similarity ratio R, such as in Eqs. (13) and (15), can be determined by using following equation:

$$R=\frac{p+{p}_{\mathrm{s}}}{\overline{p }+{p}_{\mathrm{s}}}=\frac{{p}_{\mathrm{st}}+{p}_{\mathrm{s}}}{{\overline{p} }_{\mathrm{st}}+{p}_{\mathrm{s}}}=\frac{\left(p+{p}_{\mathrm{s}}\right)(1+\frac{{\eta }^{2}}{{M}^{2}})}{{\overline{p} }_{\mathrm{st}}+{p}_{\mathrm{s}}}$$
(18)

Following authors’ previous studies on sub-loading elastoplasticity [59, 60, 62], stress ratio η for sub-loading surface and yield surface is identical, which can be served as a link of the stress point on the sub-loading surface with that on the yield surface.

3.3 Flow rule

The associated flow rule (f = g) is employed in this study to achieve plastic deviator strain as:

$$\mathrm{d}{\varepsilon }_{\mathrm{d}}^{\mathrm{tp}}=\frac{2\eta }{{M}^{2}-{\eta }^{2}}\mathrm{d}{\varepsilon }_{\mathrm{v}}^{\mathrm{tp}}$$
(19)

The increment for elastic deviator strain is written as Eq. (20):

$$\mathrm{d}{\varepsilon }_{\mathrm{d}}^{\mathrm{e}}=\frac{2(1+\nu )}{9(1-2\nu )}\frac{\kappa }{1+{e}_{0}}\frac{\mathrm{d}q}{p}$$
(20)

where ν is Poisson's ratio.

4 Model calibration

The proposed constitutive framework (i.e. the UTUH model) for unsaturated soils involves thirteen model parameters, as shown in Table 1. With the exception of eleven mechanical parameters (i.e. e(0) and ten BBM model parameters), there are two creep parameters. e(0) is the starting point of the INCL0, and it determines the position of the INCL0, which takes a role of tracing loading history together with the overconsolidation parameter R. Here, the determination of the INCL0 (or e(0)) is presented.

Table 1 Parameters for the UTUH model

The determination of the INCL0 (or e(0)) is as follows. According to the isotach concept (e.g. the loading rate lines shown in Fig. 1) proportional loads, e.g. 250, 500, 1000 and 2000 kPa, and constant duration (e.g. Δt = 1.5 days) of the adjacent loads will lead to constant deformation (i.e. Δe = −0.0832) under constant suction (i.e. s = 0 kPa), and then a straight line (e.g. the NCL0, parallel to the INCL0) with the slope of λ(0) in the semilogarithmic scale during an isotropic compression test [60]. Similarly, loading lines with different loading rates run parallel to each other under the same suction, and there must be a loading line with the slope of λ(0) passing through the initial point under a certain loading rate. Hence, the linear relationship can be determined as:

$$\frac{\mathrm{d}p}{p}=-\frac{\mathrm{d}e}{\lambda \left(0\right)}$$
(21)

With Eq. (13), total volumetric strain under isotropic condition and the saturated state can be written as:

$$\mathrm{d}{\varepsilon }_{\text{v}}=\left(\frac{\kappa }{1+{e}_{0}}+\frac{\lambda \left(0\right)-\kappa }{1+{e}_{0}}\frac{{M}^{4}}{{M}_{\mathrm{f}}^{4}}\right)\frac{\mathrm{d}p}{p}+\frac{\beta \left(0\right)}{1+{e}_{0}}\frac{{R}^{\mathrm{\alpha }}}{{t}_{0}}\mathrm{d}t$$
(22)

Hence, substituting Eq. (21) into Eq. (22) results in the relationship between de/dt and R under the saturated state as:

$$ \frac{{{\text{d}}e}}{{{\text{d}}t}} = \beta (0)\frac{{\lambda \left( 0 \right)}}{{\kappa - \lambda \left( 0 \right)}}\frac{{R^{\alpha } }}{{t_{0} }}\left( {1 - \frac{{M^{4} }}{{M_{{\text{f}}}^{4} }}} \right)^{{ - 1}} $$
(23)

By Eq. (23), the initial value of R can be determined (i.e. 0.85) with Δe = −0.0832 and Δt = 1.5 days. According to the definition of R, the swelling line and an arbitrary data point (e.g. B (2000 kPa, 1.0979) in Fig. 5), a point on the INCL0 (e.g. A (2352.9 kPa, 1.0977) in Fig. 5) can be obtained. Finally, with the point A and the slope λ(0), the INCL0 is determined as \(e=2.029-0.12\mathrm{ln}\frac{p}{{p}^{\mathrm{c}}}\) (i.e. e(0) = 2.029).

Fig. 5
figure 5

Determination of the INCL0

To determine the parameter b, nonlinear regression analysis can be carried out, using test date from Qin et al. [41], Ye et al. [64] and Esfandiari et al. [15]. Qin et al. [41], Ye et al. [64] and Esfandiari et al. [15] have experimentally investigated time-dependent deformation behaviour of unsaturated soils. For example, the compression index and the secondary compression index at different matric suction values were measured by using oedometer. Here, the relationships between the compression index, the secondary compression index and suction are analysed to determine the parameter b with \(\varphi \left(s\right)=\varphi \left(0\right)-b\mathrm{ln}\left(\frac{s+{p}_{at}}{{p}_{at}}\right)\), where φ(0) is ratio of coefficient of the secondary consolidation to compression coefficient when suction is 0 (saturated soils). Finally, b is identified, as shown in Fig. 6.

Fig. 6
figure 6

Comparison between data [15, 41, 64] and model calibration

5 Model performance

In this section, the performance of the proposed UTUH model is illustrated with four numerical examples. In the first example, the performance of hypothetical unsaturated specimens with constant suction and different loading rates is investigated under isotropic conditions. The second example is used to assess the effect of stress history under isotropic conditions, and it includes isotropic loading tests with different loading rates for unsaturated soil samples to achieve different stress history. The third example and the fourth example are employed for investigating response of unsaturated soils to different wetting rates under isotropic conditions and non-isotropic conditions, respectively. In both examples, different wetting rates are designed when unsaturated soil samples are subjected to wetting. Table 2 shows the parameter values used in model performance.

Table 2 The parameter values used in model performance

5.1 Evolution of constitutive framework for unsaturated soils under isotropic conditions and constant suction

Performance of the UTUH model under isotropic conditions is studied here. The hypothetical unsaturated specimens with constant suction 0 kPa and 100 kPa are loaded under different loading rates, and the p versus e results are plotted in Fig. 7. Figure 7 shows all three characteristics, as mentioned in Fig. 1 in Sect. 2.1. Besides, it can also be found that starting at the same initial point under the same suction, the curves corresponding to faster volume strain rates gradually merge into the INCLs with the elapsing time of loading, which also means that they are located at higher position. However, for those with slower volume strain rates, they show totally reverse trends. The slower volume strain rates they have, the farther away from the INCLs they are, and the higher overconsolidation degree they have. With a certain loading rate (i.e. 1.1 × 10–3), a linear loading line with the slope of λ(0) passing through the initial point can be determined under constant suction 0 kPa in the semilogarithmic scale. Figure 7 also presents that the proposed UTUH model can describe relaxation.

Fig. 7
figure 7

Evolution of the UTUH model under isotropic conditions and constant suction

For unsaturated soils (e.g. s = 100 kPa), if 0.55 is taken as the target value of void ratio (see Fig. 7), it is found that to reach it immediately (rate is infinity) from the initial point (e = 0.72), the required pressure is higher than 4600 kPa. However, net mean stress is less than 3500 kPa, when following the rate line 4.4 × 10–5 in as much as time induced significantly void ratio reduction. When we consider constant net mean stress (e.g. p = 3000 kPa), there is significant void ratio compression for the rate line 4.4 × 10–5 compared to that for the infinite rate line owing to the contribution from time. Figure 7 represents deformation of unsaturated soils, which experience different loading time and different loads under constant suction.

5.2 Behaviour of unsaturated soils with different stress history under isotropic conditions

In the case II, three scenarios are chosen for investigating the effect of stress history on the performance of the UTUH model under the isotropic state. For example, the first scenario involves the drying path, the loading path and the wetting path, as shown in Fig. 8a, and all the three paths are time-independent. The drying path starts under the saturated state and at confining stress 100 kPa, and suction increases from 0 to 200 kPa during drying. Then the hypothetical soil sample is isotropically compressed at constant suction (s = 200 kPa) from 100 to 6 × 104 kPa. Finally, suction reduction begins from 200 kPa under controlled net mean stress 6 × 104 kPa until to 0 kPa. The first scenario is taken as the benchmark, and the whole process is time-independent while the drying paths for the second scenario and the third scenario, and the loading path for the third scenario involves the time effect (i.e. the drying rate is 1 × 10–2 kPa/min and the loading rate is 1 × 10–6/min), as shown in Fig. 8a. The prediction results for the three scenarios are given in Fig. 8. For instance, Fig. 8a shows stress paths, and void ratio-net mean stress relationships are plotted in Fig. 8b. Similarly, relationships between void ratio and time and suction are shown in Fig. 8d, e, respectively. Figure 8c displays the changes of void ratio (i.e. the gaps between the INCL0 and loading lines in Fig. 8b) with net mean stress.

Fig. 8
figure 8

Prediction results for unsaturated soils with different stress history under isotropic conditions

Starting at a common initial point, three hypothetical soil samples are dried from the saturated condition to the unsaturated state, and suction increases from 0 to 200 kPa. During the drying process, there is only slight elastic deformation in the first scenario since the increase of suction is not outside the elastic region. However, the two hypothetical soil samples for the second scenario and the third scenario not only experience the increase of elastic strain, but also viscous deformation (see Fig. 8e). It reaches the peak at the end of drying (see the gap in Fig. 8c, e). Figure 8d shows that soils are prone to deformation at the saturated state, and reduction of the effect of time on deformation can also be observed given that slopes for dashed lines in Fig. 8d descend as suction increases. This is because drying is favourable to the increase of soil stiffness, and then unsaturated soils build up resistance to viscous deformation. After the drying process, the gap between the solid line and dashed lines induced by the creep becomes increasingly narrow as p increases to a high value, as shown in Fig. 8c. During the loading, Fig. 8b shows that the compression curves for the first scenario and the second scenario merge into the INCL200 while that for the third scenario approaches to the rate line (i.e. 1 × 10–6/min). Hence, there is a gap between them, and the combination of p-loading and the t-increasing results in the largest gap at the end of the loading process (see Fig. 8b–d). During wetting, the collapses of three hypothetical soil samples begin once suction reduction starts. Finally, the wetting curve for the first scenario merges into the INCL0 (see Fig. 8b, c). However, that for the second scenario and the third scenario is below the INCL0, and the third hypothetical soil sample experiences the largest changes in void ratio of the three due to the time effect when the wetting process is completed. Figures 8c, d show that the differences of void ratio between the second scenario and the third scenario are increasingly narrow as wetting in as much as when isotropic compression is completed, the hypothetical unsaturated soil sample has high overconsolidation degree and then poor collapsibility, which is less likely to experience deformation. This phenomenon is consistent with the case in geotechnical engineering when earth structures under the unsaturated condition undergo different loading and creep histories due to complex environments.

5.3 Behaviour of unsaturated soils with different wetting rates under the isotropic state

In the case III, the performance of the UTUH model is assessed upon the isotropic state, and response of unsaturated soils to different wetting rates (i.e. the wetting rates are negative infinity, −2 × 10–3 kPa/min and −5 × 10–4 kPa/min) is highlighted in this section. In the case III, three scenarios are adopted, and each of them includes the drying path, the loading path and the wetting path (see Fig. 9a). During drying, the three hypothetical saturated samples for the three scenarios are dried with the constant drying rates 2 × 10–3 kPa/min at confining stress 100 kPa until suction climbs to 200 kPa. Then increasing net mean stress from 100 kPa to 6 × 104 kPa is applied to three hypothetical unsaturated samples with constant suction 200 kPa under the isotropic state, and the loading rate is 1 × 10–6/min. Once loading is completed, the wetting process stars under controlled net mean stress 6 × 104 kPa, and suction decreases from 200 to 0 kPa with the three wetting rates, negative infinity, −2 × 10–3 kPa/min and −5 × 10–4 kPa/min corresponding to the first scenario, the second scenario and the third scenario, respectively, as shown in Fig. 9a. It should be noted that except for the wetting path for the first scenario, all paths for the three scenarios are time-dependent. The simulation results of the three scenarios are plotted in Fig. 9. For example, Fig. 9a gives stress paths, and void ratio–net mean stress relationships are shown in Fig. 9b. Besides, void ratio vs. time, void ratio vs. suction and the overconsolidation parameter R vs. suction are plotted in Fig. 9d–f, respectively. Figure 9c displays the changes of void ratio with suction.

Fig. 9
figure 9

Prediction results for unsaturated soils in different wetting rates under isotropic conditions

Beginning at the initial stress point, void ratio of the three hypothetical soil samples for the three scenarios shows a downward trajectory (see Fig. 9b–e) during drying due to slight elastic deformation induced by the increase of suction, and continuous time-dependent deformation. The decrease of R from the initial value 1 to a very small value (see Fig. 9f) implies the increasing overconsolidation degree. Drying leads to the gradual increase of unsaturated soil stiffness. This increase and poor compressibility for unsaturated soil owing to the increasing overconsolidation degree cause descending slopes of creep, as shown in Fig. 9d. During time-dependent isotropic compression, all the three compression curves gradually merge into the rate line (i.e. 1 × 10–6/min) with the same slopes of INCL200, as shown in Fig. 9b. The gap between the two lines (INCL200 and 1 × 10–6/min) indicates the time effect. Before time-dependent isotropic compression is completed, the trajectories of void ratio, net mean stress, time, suction and the overconsolidation parameter are the same (see Fig. 9a–f) under the three paths for the three scenarios because of the same loading strategy. Nevertheless, once wetting stars, void ratio, time and the overconsolidation parameter follow a unique path now that there are three wetting rates (i.e. the wetting rates are negative infinity, −2 × 10–3 kPa/min and −5 × 10–4 kPa/min). During wetting steps, the collapses and creep of three hypothetical soil samples induce continuous volumetric deformation. Figure 9b–e shows that when wetting is completed, the third hypothetical soil sample for the third scenario experiences the largest void ratio reduction and the smallest overconsolidation parameter of the three scenarios due to the lowest wetting rates and then the longest creep time of three. Figure 9d indicates that the compression curves corresponding to the second scenario and the third scenario own slightly increasing slopes as wetting since wetting results in reduction of unsaturated soil stiffness and the increase of the viscous parameter, which pose an adverse effect on structures built on unsaturated soils in practice.

After time-dependent isotropic compression, void ratio of unsaturated soils is 1.28, and unsaturated soils are overconsolidated (see Fig. 9f). If 0.8 is taken as the final value of void ratio, Fig. 9d shows that to reach the reduction in void ratio (Δe = 1.28–0.8 = 0.48) immediately (t = 0), the required suction changes are 183 kPa (Δs = 200–17 = 183 kPa) while suction reduction for the second scenario and the third scenario are 176 kPa (Δs = 200–24 = 176 kPa) and 168.5 kPa (Δs = 200–31.5 = 168.5 kPa), respectively, for achieving the same void ratio changes, as shown in Fig. 9e. This is because creep in 61 days (the second scenario) and 234 days (the third scenario) makes contributions to void ratio reduction (see Fig. 9d). Consider other case. Compared with void ratio compression for the third scenario (Δe = 1.28–1.16 = 0.12), that for the second scenario reaches to a very high value (Δe = 1.28–0.81 = 0.47) on conditions that the final creep time is kept as unchanged (e.g. 100 days, see Fig. 9d) because unsaturated soils for the second scenario experience the more serious collapses due to suction reduction (Δs = 200–28.5 = 171.5 kPa, see Fig. 9e). Hence, it can be seen that Fig. 9 represents deformation of unsaturated soils subjected to different creep time and different wetting loads. Understanding and predicting deformation of unsaturated soils subjected to different creep time and different wetting loads have highly practical significance in geotechnical engineering since earth structures under the unsaturated condition inevitably encounter different natural or man-induced saturation processes, which take different periods for changeable weather (e.g. quick wetting due to a heavy rainstorm and slow wetting due to long-last small rain).

5.4 Behaviour of unsaturated soils with different wetting rates under the triaxial stress state

Consider the triaxial stress state. Three scenarios are employed here for illustrating the application of the UTUH model. Stress paths for the three scenarios successively involve shearing steps, wetting steps and shearing steps. Unsaturated soil samples are subjected to shearing under constant suction 100 kPa until q/p reaches to 1.26. There followed wetting steps, suction reduces from 100 to 0 kPa, and wetting rates for the three scenarios are negative infinity, −2 × 10–2 kPa/min and −1 × 10–3 kPa/min, respectively. The loading process ends with shearing steps. Figure 10 shows prediction results for the three scenarios. Figure 10a presents stress and strain response, and Fig. 10b gives deviator strain–time relationships. Besides, volumetric strain vs. time is shown in Fig. 10c.

Fig. 10
figure 10

Prediction results for unsaturated soils in different wetting rates under the triaxial stress state

During isotropic compression, time-dependent wetting can cause only volumetric strain. However, when unsaturated soils are under the triaxial stress state, it may not only induce volumetric strain, but also deviator strain under constant deviator stress due to the associated flow rule, as shown in Fig. 10. Figure 10 also presents that the slower wetting rates indicate the larger deviator strain and the larger volumetric strain for unsaturated soils since the slower wetting rates mean the greater deformation contribution from time. Figure 10a suggests that overconsolidation behaviour of unsaturated soils (i.e. softening of strain and dilatancy of volume) becomes increasingly apparent, and the peak for q/p increases as wetting rates decrease because the slower wetting rates there are, the greater overconsolidation degree the unsaturated soil samples have. Figure 10b shows that when the first shearing process is completed, deviator strain of unsaturated soil samples is 0.024. To reach the target value 0.08 of deviator strain, the third scenario needs the longer creep time as suction reduction in the third scenario makes smaller contributions to deviator strain. On the other hand, if creep time is kept as constant (e.g. 3 days, see Fig. 10b), there is a higher value of deviator strain of the unsaturated soil sample for the second scenario as it experiences the more serious collapses due to suction reduction. Figure 10c shows the similar behaviour for volumetric strain.

6 Model validations

Qin et al. [41], Banerjee and Puppala [7], Banerjee [6] and Patil et al. [36] have taken experimental effort to examine time-dependent behaviour of unsaturated soils. In this study, their test results are utilized to evaluate validity of the developed UTUH model through four cases. Corresponding parameter values are listed in Table 3.

Table 3 The parameter values used for validations

In case I, oedometer test results from Qin et al. [41] are used for validations. In the test, Qin et al. [41] have investigated strain rate-dependent behaviour of unsaturated GMZ01 bentonite under controlled suctions and CRS conditions. For example, the first group consists of one sample with constant suction (i.e. s = 0 MPa) and initial void ratio 1.105, and it experiences compression from initial vertical stress 0.1 MPa to final vertical stress 20 MPa under CRS 1 × 10–5, 1 × 10–6, 1 × 10–7 and 1 × 10–5/s, sequentially. The experimental date and the model prediction outcomes are marked as the green points and the green solid line in Fig. 11, respectively. The prediction results agree well with the test date. Owing to the different strain rates, there is step variation of the loading curves, and lower strain rates lead to smaller void ratio in Fig. 11. The second group consists of three samples with controlled suction (i.e. s = 9 MPa) and initial void ratio 0.745, and they are subjected to compression from initial vertical stress 0.08 MPa to final vertical stress 35 MPa under single-stage CRS (i.e. 1 × 10–5 and 1 × 10–6/s) and stepwise CRS (i.e. 1 × 10–5, 1 × 10–6, 1 × 10–7, 1 × 10–5, 1 × 10–6 and 1 × 10–7/s), respectively. The comparison between test date and prediction results including the single-stage CRS and the stepwise CRS is also shown in Fig. 11. A higher level of suction (i.e. s = 110 MPa) is applied in the third group consisting of one sample with initial void ratio 0.564. It is compressed from initial vertical stress 0.1 MPa to final vertical stress 55 MPa under stepwise CRS (i.e. 1 × 10–6, 1 × 10–7 and 1 × 10–5/s). The test date and the model prediction results are represented as the magenta points and the magenta solid line in Fig. 11, respectively. The proposed UTUH model captures step variation of the loading curves due to the different strain rates. Figure 11 also indicates that with the same suction (s), loading lines with different strain rates run nearly parallel to each other, the slope for loading lines begins its ascent with the decrease of suction, and the gap between adjacent loading lines becomes increasingly narrow with the increase of suction, which means the viscous parameter β increases during wetting. Figure 11 shows that under a controlled suction (e.g. s = 0, 9 or 110 MPa) and a constant vertical stress, void ratio gradually decreases as strain rates steadily drop from 1 × 10–5/s to 1 × 10–7/s, indicating that the increase of loading time causes additional deformation of GMZ01 bentonite.

Fig. 11
figure 11

Validation case I: CRS test results for soils under controlled suctions (data after Qin et al. [41])

Banerjee and Puppala [7] systematically investigated the effect of loading rates on the strength and deformation of unsaturated soils via the conventional triaxial compression (CTC) tests under two suction levels, and their tests are represented by case II in this paper. In their CTC experiment, the first series consists of four samples, and the four samples with the consistent suction level (i.e. s = 0 kPa) are subjected to shearing in initial confining pressures 400 kPa under four rates (i.e. 5 × 10–2, 1 × 10–1, 2.5 × 10–1 and 5 × 10–1%/min) until axis strain reaches to 0.2. Figure 12 shows corresponding comparison results between test date and model prediction. The proposed UTUH model well captures stress and strain features in the measured behaviour although it overestimates the peak strength of soils. Figure 12 suggests that under a constant suction level (i.e. s = 0 kPa), the decrease of loading rates causes the decline of the peak shear strength and the increase of volumetric deformation of soils due to the gradually reinforced time effect. For example, the peak shear strength for 5 × 10–1%/min is 956 kPa while that for 5 × 10–2 is 838 kPa.

Fig. 12
figure 12

Validation case II: CTC test results for soils with initial confining pressures 400 kPa and constant suction (i.e. 0 kPa), and data after Banerjee and Puppala [7]

In CTC experiment conducted by Banerjee and Puppala [7], the second series consists of three samples with the same suction (i.e. s = 250 kPa), and they experience shearing in initial confining pressures 400 kPa under three rates (i.e. 1 × 10–3, 3 × 10–3 and 1 × 10–2%/min) until axis strain reaches to 0.12. Figure 13 gives corresponding test date and model prediction outcomes. The UTUH model nicely predicts stress and strain features in the measured behaviour. Figure 13 also shows that under constant suction (i.e. s = 250 kPa), the time effect results in the decrease of the shear strength and the increase of volumetric deformation of unsaturated soils. For example, when strain rates decrease from 1 × 10–2 to 1 × 10–3%/min, there is also the decline of the peak shear strength from 1376 to 1155 kPa.

Fig. 13
figure 13

Validation case II: CTC test results for soils with initial confining pressures 400 kPa and constant suction (i.e. 250 kPa), and data after Banerjee and Puppala [7]

Banerjee [6] studied responses of strength and deformation behaviour of unsaturated soils to different loading rates via CTC experiment under two suction levels (i.e. s = 0 and 250 kPa), and it is represented by case III in this section. In this CTC experiment, the first series includes four samples, and the four samples with the same suction (i.e. s = 0 kPa) are sheared in initial confining pressures 400 kPa under four rates (i.e. 1 × 10–2, 5 × 10–2, 2.5 × 10–1 and 5 × 10–1%/min) until axis strain reaches to 0.2. Figure 14 presents corresponding comparison results between test date and model prediction. The UTUH model well predicts stress and strain features in the observed behaviour. Figure 14 shows that under a constant suction level (i.e. s = 0 kPa), the increase of loading rates causes the increase of the peak shear strength and the decrease of volumetric deformation of soils due to the gradually weakened time effect. For instance, the peak shear strength for the sample with rates 1 × 10–2%/min is 990 kPa while that for the sample with rates 5 × 10–1%/min is 1153 kPa.

Fig. 14
figure 14

Validation case III: CTC test results for soils with initial confining pressures 400 kPa and constant suction (i.e. 0 kPa), and data after Banerjee [6]

In CTC experiment performed by Banerjee [6], the second series includes three samples, and the three samples with the same suction (i.e. s = 250 kPa) experience shearing in initial confining pressures 400 kPa under three rates (i.e. 1 × 10–3, 3 × 10–3 and 1 × 10–2%/min) until axis strain reaches to 0.16. Figure 15 gives corresponding comparison results between test date and model prediction. The UTUH model gives nice prediction of stress features in the measured behaviour. It is found in Fig. 15 that under the same suction (i.e. s = 250 kPa), the time effect is unfavourable to the increase of the peak shear strength of unsaturated soils, but causes the increase of volumetric deformation. For example, the peak shear strength for the sample with rates 1 × 10–2%/min is 1385 kPa while that for the sample with rates 1 × 10–3%/min is 1215 kPa. Besides, Figs. 14 and 15 present that soils with a constant strain rate (i.e. 1 × 10–2) experience the strength increase with the increase of suction. When suction increases from 0 to 250 kPa, the peak shear strength increases from 990 to 1385 kPa under the same strain rate (i.e. 1 × 10–2). This is because drying is favourable to the increase of soil stiffness, and then soils build up resistance to shearing.

Fig. 15
figure 15

Validation case III: CTC test results for soils with initial confining pressures 400 kPa and constant suction (i.e. 250 kPa), and data after Banerjee [6]

Patil et al. [36] studied strain–stress behaviour of unsaturated soils with different loading rates. In their CTC experiment, three samples with the consistent suction level (i.e. s = 500 kPa) are subjected to shearing in initial confining pressures 300 kPa under three rates (i.e. 2.9 × 10–3, 8.6 × 10–3 and 1.4 × 10–2%/min). The test date is plotted together with simulated outcomes in Fig. 16. The UTUH model gives a satisfactory description of stress–strain behaviour of unsaturated soils. Figure 16 shows that under a constant suction level (i.e. s = 500 kPa), the decrease of loading rates causes the decline of the peak shear strength and the increase of volumetric deformation of unsaturated soils. Besides, during shearing, unsaturated soils show hardening behaviour as deviator stress increase until they reach the peak points. Once the peak points are reached, deviator stress for unsaturated soils decreases, and finally it stabilizes at constants, which is known as softening behaviour.

Fig. 16
figure 16

Validation case IV: CTC test results for soils with initial confining pressures 300 kPa and constant suction (i.e. 500 kPa), and data after Patil et al. [36]

In this section, it is found that comparison in the four cases is satisfactorily in agreement with the test results. The good agreements indicate that the proposed UTUH model is capable of reproducing time-dependent aspects of unsaturated soils, e.g. the effect of loading rates on the shear strength and volumetric deformation, and the strength increases due to the increase of suction. Furthermore, the developed UTUH model can describe three key characteristics, as mentioned in Fig. 1 in Sect. 2.1.

7 Conclusions

In this study, INCLs for unsaturated soils is introduced to determine creep time and overconsolidation states of unsaturated soils. With the INCLs, the increment of viscous deformation in an equivalent form of the overconsolidation parameter is derived. Subsequently, an isotropic EVP constitutive framework for unsaturated soils is proposed by combining viscous deformation with mechanical and hydraulic deformation through overconsolidation parameter.

A yield surface and a sub-loading surface are introduced, and an extension to a triaxial stress state is built in the space of mean effective stress (p), suction (s), deviator stress (q) and time variable (\(\widetilde{t}\)). Then, the way to determine model parameters and performance of the proposed model are demonstrated. The proposed EVP constitutive model, i.e. the UTUH model, can characterize the mechanical, time-dependent and hydraulic behaviour of unsaturated soils, including the joint effects of mechanical/hydraulic loading and creep histories. The response of unsaturated soils to different wetting rates is well assessed with the UTUH model under isotropic conditions and non-isotropic conditions. Experimental results from literature have been employed to validate the UTUH and confirmed that the UTUH model is capable of reproducing time-dependent aspects of unsaturated soils, e.g. the effect of loading rates on the shear strength and volumetric deformation, and the strength increases due to the increase of suction.

In the future work, the UTUH model can be implemented in finite element code to solve boundary value problems by using existing FE implementation strategies [21] to assess the performance of earth structures build on/in unsaturated soils in future work. In addition, it should be acknowledged that the behaviour of unsaturated soils could be very complex owing to their composition and stress conditions. More investigation concerning hydro-mechanical coupling and anisotropic behaviour as well as brittle response of unsaturated soils with high suction can be performed in a future study through incorporating more advanced models, such as hydro-mechanical models [28], anisotropic models [63] and models feature fracture mechanics [37].