Disturbed state concept and non-isothermal shear strength model for unsaturated soils

Shear strength of unsaturated soils is an important engineering property that is required for addressing geotechnical problems, the prediction of which remains to be a challenging task for design engineers due to the complex interaction problem. This study presents a new shear strength equation based on the micromechanical model and the disturbed state concept for unsaturated soils. The original point of this study is considering the solid contact area ratio which was neglected in most of the existing equations. Using the proposed model, the non-linear relationship between the matric suction, saturation degree, and the shear strength of unsaturated soils are described. Validation of the shear strength model was verified against the experimental data and several current models on six different types of soils. The results indicate that the proposed model has a good performance in predicting the shear strength of unsaturated soils, and generally is better than other existing models. In response to varying climatic conditions, the analytical model was then extended to consider the effect of temperature on the shear strength of unsaturated soils. The comparison between predicted and measured results was carried out on compacted silt for three different temperatures. The results show that the proposed model is capable of accurately predicting changes in unsaturated shear strength as a function of temperature.


Introduction
Geotechnical engineering research in general and unsaturated soil mechanics remain a crucial part of climate change mitigation and adaptation. The temperature increase due to climate change has driven the evaporation as well as changes in subsurface conditions. In these situations, the water amount in soils is constantly changed and soils tend to change from saturated to unsaturated state. As a result, a significant portion of the soil is involved in unsaturated conditions, particularly in the surface soil layers, seasonal areas, soils above the water table, and compacted soils. On the other hand, it is also well-known that the stability of geotechnical structures depends strongly on the accuracy in predicting the shear strength of soils. In this context, the understanding of shear strength principles of unsaturated soils becomes a topic of interest in addressing geotechnical problems. It is, therefore, crucial to quantify the shear strength of unsaturated soils as well as to be able to respond to the effect of climate change.
The shear strength of unsaturated soils is a highly challenging topic mainly due to their complex multi-phase interaction nature. In practice, the unsaturated shear strength is often predicted by using the two independent stress state variables, named as net normal stress (σ − u a ) and matric suction (u a − u w ). To have a better understanding of unsaturated shear strength, several researchers have efforted to conduct experimental campaigns over the last decades using different test types such as the triaxial test (Blatz and Graham 2003;Houston et al. 2008;Goh et al. 2014;Kim et al. 2018;Satyanaga and Rahardjo 2019), modified triaxial test (Estabragh and Javadi 2008;Patil et al. 2016a, b;Zheng et al. 2020), direct shear test (Hamid and Miller 2009;Nam et al. 2011;Schnellmann et al. 2013;Khaboushan et al. 2018). The experimental results emphasized that the unsaturated shear strength of the soil increases with increasing net normal stress and matric suction, following a nonlinear relationship, yet the level of the nonlinearity depends on the applied net normal stress. The nonlinearity of the shear strength envelope is related to the changes in the degree of saturation of the soil as soil suction changes. Finally, the critical state was also observed and confirmed for unsaturated soils. The above observations should be considered in the development of any analytical models (Pham and Sutman 2021).
Unfortunately, the experimental determination of the shear strength of unsaturated soils is often expensive, timeconsuming, and complex. The challenges related to experimental approaches render the numerical methods as a prominent tool in the framework of unsaturated soils. Being a multiphase system, the response of unsaturated soils is based on the grain-scale micromechanical mechanism and therefore shows a particular behaviour compared to other types of materials (Lloret-Cabot et al. 2017). As a result, several different numerical methods were used to study the shear behaviour of unsaturated soils, such as the discrete element method (Jiang et al. 2004;Cai et al. 2020;Kim and Park 2020), finite difference method (Gong et al. 2018), and finite element method (Cho and Lee 2001;Abed and Vermeer 2009). However, the numerical studies related to the shear behaviour of unsaturated soils are still limited due to the complexity of the problem.
In this context, the theoretical framework with mathematical equations also becomes an important tool to predict the shear strength of unsaturated soils. This is particularly true at the preliminary design stage of a project. It is worthy that the general principles associated with the shear strength of saturated soils are often used to apply to unsaturated soils. For engineering purposes, several different equations were proposed to predict the shear strength of unsaturated soils (Bishop et al. 1960;Fredlund et al. 1978;Toll 1990;Vanapalli et al. 1996;Oberg and Sallfors 1997;Khalili and Khabbaz 1998;Miao et al. 2002;Lee et al. 2005;Matsushi and Matsukura 2006;Kayadelen et al. 2007;Zhai et al. 2019). In most cases, the estimated equations are based on the saturated shear strength parameters and the soil-water characteristic curve (SWCC) of the unsaturated soil. Among the existing shear strength equations, most of them are empirical and phenomenological by using fitting parameters. However, the empirical equations are limited to a narrow range of soils employed in the corresponding experiments. Another limitation is that these equations usually require a significant number of input parameters that are not directly available from experimental studies, or they are time-consuming to be identified by laboratory tests. In practice, these methods may produce different results for the same soil sample (Vanapalli and Fredlund 2000;Ye et al. 2010;Patil et al. 2017;Tang et al. 2019;Pham and Sutman 2021;Pham 2022).
It is therefore always desirable to develop a shear strength equation based on the micromechanical theory considering a multi-phase interaction system for unsaturated soils. In this study, the micromechanical interaction model combined with disturbed state theory is first presented to provide a new equation for predicting the shear strength of unsaturated soils. The proposed model is compared with several existing equations and the experimental results using six published data sets to investigate its validation in the next section. In response to varying climatic conditions, the proposed model is extended to apply to a non-isothermal condition that allows estimating the unsaturated shear strength with temperature changes.

Micromechanical analysis of unsaturated soils
In this section, the micromechanical equilibrium model presented in Pham and Sutman (2021) is refined by the inclusion of a solid contact area ratio by considering the disturbed state concept. In the development of the shear strength model, the following principles are applied: a) An unsaturated soil can be considered a mixture of three phases, which include solids, pore water, and pore air. b) The shear strength of unsaturated soils is the result of the interaction between the three phases, and the breaking of the bonds. c) The models of the individual phases are assembled, using micromechanical equilibrium, to predict the overall shear strength of unsaturated soils.
Therefore, the following analysis starts with the micromechanical stress equilibrium model of three phases to develop a new effective stress equation for unsaturated soils.

Unsaturated effective stress considering the disturbed state
When soil particles are loaded, they tend to reorient and come into contact with one another. The actual example is simplified by postulating a soil made up of homogeneous rigid spheres to better explain the interaction between phases. However, it should be noted that the solution of the analytical model in this paper is established based on the volume relationship and force equilibrium on the total section area without being associated with particle size or particle arrangement. Figure 1a shows the contact of two soil grains in an unsaturated soil. The area of contact between two soil particles is noted as A c , and the area of the two water membranes embracing the contact area is A m w1 and A m w2 . The load transfer between the two grains is obtained partly through the inter-grain contact area A c , and partly through the menisci water areas, which can be written as: Most natural soils, on the other hand, have a bimodal structure with a microstructure and a macrostructure. The majority of the macrostructure and the entire microstructure stay saturated when the matric suction range applied to the soil is low. Meanwhile, as suction increases, the saturated soil zone decreases, with only the microstructure saturated (Romero et al. 2001;Rojas 2008;Alonso et al. 2010). As a result, the unsaturated soil is divided into two parts: a saturated portion with "free water" around the soil particles, and an unsaturated fraction with solid particles joined together by water menisci. Because there are certain contact areas between soil particles in free water, it can bear an external load together with the soil skeleton (Karube and Kawai 2001). A water distribution function f is introduced to consider the effect of free water, which was previously overlooked by existing models, and the total water area in an unsaturated soil sample must therefore be: The interaction force system in unsaturated soils is depicted in Fig. 1b. Within the framework of multi-phase continuum mechanics, the equilibrium between internal and external forces for each phase can be described as follows: where P is the total external load applied on the unsaturated soil, P c is the force transferred through the grain contact, P w is the total force carried by pore water, P a is the force carried by pore air. Figure 1c depicts the stress distribution in an unsaturated soil cross-section. It is worth noting that the solid particles make contact with all three phases; therefore, there are two menisci on either side of each contact point. As a result, the force carried by pore water in which the liquid phase reacts, pore air in which the air phase reacts, and solid contact in which the solid phase reacts may be represented as: where P m w1 and P m w2 are the forces carried by water in two menisci surrounding each contact point respectively, u w is the pore water pressure, u a is the pore air pressure, A is the 1 Stress equilibrium analysis model: a representative contact area of unsaturated soils, b force interaction analysis, c stress distribution in unsaturated soil total section area, A c is the solid contact area between the two soil grains, A w is the total water area, f is a factor to consider the combined effect of free water and menisci water in unsaturated soils, is total stress, c is solid contact stress, ′ is effective stress that represents the equivalent contact stress over the total section area.
In soil mechanics, the term "effective stress" has been adopted in lieu of solid contact stress when calculating the shear strength of the material. This is because soils are a mixture of continuous and discontinuous parts, and the failure state of soils depends on the average stress distribution over the material elements (nonlocal), not at a point (local).
Substituting Eqs. (4) to (7) back into Eq. (3) gives: Dividing both sides of Eq. (8) to the total cross-section area (A) results in: Rearranging Eq. (9) gives where D is the disturbance function, defined as the ratio of solid contact area to total area: The effective stress equation for unsaturated soils is given by rearranging Eq. (10) as follows: It should be noted that the unsaturated pore water pressure can be calculated by multiplying the pore water pressure at the saturated state with the saturation degree of soils. Rojas (2008) suggested an equivalent stress model for bimodal structured soil (including micro and macro voids), in which an unsaturated soil is assumed to be divided into two fractions: a saturated fraction, in which soil particles are surrounded by water, and an unsaturated fraction, in which solid particles are linked together by water menisci. Here, a similar technique is utilized, with free water representing the saturated fraction and menisci water representing the unsaturated fraction. The following is the water distribution

Water distribution function in unsaturated soils
function that describes the volume relationship between free water, menisci water, and global water: where V f w is free water volume, V m w is menisci water volume, and V s is solid volume.
It should be noted that in the case where soils approach a saturated state, V f w + V m w = V w , and f = V/V w = 1/θ. Otherwise, in the case where soils approach very high suction, V m w = V w , and f = 1. If it is admitted that the value of suction is the same everywhere in the sample, then it can be concluded that all saturated zones are surrounded by menisci of water showing the same radius of curvature as the unsaturated zones. Considering the effective degree of saturation, Eq. (13) can be expressed as: where e is normalized volumetric water content, S e is effective saturation degree.
On the other side, it should be noted that the matric suction is a derivative of the Helmholtz free energy, in which the water distribution function depends on the saturation degree of soils (Lamborn 1986;Aubeny and Lytton 2003;Mun 2005). Equation (14) is in good agreement with the curve of Helmholtz's free energy theory, as illustrated in Fig. 2. The top and lower boundaries of the water distribution function in unsaturated soils are provided by Eq. (14). f = 1/ e for fully saturated soils (zero suction) while f = 1 for dry soils (very high soil suction).

Determination of ratio A m w ∕A
If pores are randomly distributed in a homogeneous isotropic material, the areas of water, air, and solids appearing in a cross-section area can be converted to the volumetric distribution of the phases in the following form: where n is porosity of soils, where S r = residual degree of saturation.

Disturbance function
The disturbance is the term used to describe the deviation of the current deforming state concerning the initial state (reference state) of the material (Desai and Wang 2003). The disturbance, thus, also represents the variation in the density of soils. Assume that initially, the void volume is higher and that under compressive loading it decreases due to compaction. The change in the void volume of soil can be proportional to the solid-to-solid contact volume, which increases during deformation. The total solid contact volume for a given stage is derived by multiplying the contact area between two solids by the number of interaction points (Fig. 3). The total solid contact volume, therefore, can be expressed as follows: where V c is total solid contact volume, V ci is the contact volume between two solids, is the contact angle, R is the particle radius, t is the contact thickness, and n s is the number of contact points. The definition of contact angle is illustrated in Fig. 1a.
On the other side, the following formula serves to convert between two-dimensional and three-dimensional relationships: The disturbance function thus can be defined by replacing Eq. (18) in Eq. (19) as follows: where V si is the volume of a solid particle ( V si = 4 R 3 ∕3 ), e is the void ratio of the soil sample.
In this study, the contact angle at the loosest state is assumed to be approximately 10° for sands and 25° for clays, and contact thickness is considered t = 2R. Therefore, the approximate solution for the disturbance function is implied as follows: On the other hand, if the result of the shear test at a saturated state is available, the disturbance function can be derived empirically. It is possible to calculate the empirical disturbance function as follows: where u w is pore-water pressure, m is measured shear stress at a saturated state, c′ is the effective cohesion, ϕ′ is the friction angle of soils.

Discussion of Terzaghi effective stress equation
Substituting Eqs. (14) and (15) back into Eq. (12) give a general form of effective stress as follows: For unsaturated soils, Eq. (24) represents a general form of effective stress. Considering two special cases where soils are completely saturated, and soils are completely dry: If soils are completely saturated ( u a = u w ), Eq. (24) can be revised as: If soils are completely dry ( A m w = 0 ), Eq. (24) can be rewritten as: On the other side, the Terzaghi effective stress equation is also well-known as follows: By comparing Eqs. (25) and (27), it can be observed that the main difference between the proposed equation and Terzaghi's effective stress equation for saturated soils is the term disturbance (D). If the solid contact area is neglected, Eq. (25) becomes an identical one of Terzaghi's effective stress equation. Unlike saturated soils, however, the pore water pressure in unsaturated soils causes both local and non-local actions (Jennings and Burland 1962), which are reflected well in the proposed equation.

Unsaturated shear strength equation
The shear strength constitutive relationship provides a mathematical equation relating the normal and shear components of the stress tensors. To predict the shear strength, the Mohr-Coulomb failure criterion was extended to embrace unsaturated soils by Fredlund et al. (1978). The proposed shear strength equation for unsaturated soils can be written as follows: where c′ is the effective cohesion, φ′ is the friction angle of soils under saturated conditions.
Considering the proposed equation, the soil property function (SPF), which defines the relationship between shear strength and soil suction, is therefore derived as follows: Alternatively, several researchers preferred to use the term "cohesion intercept" to emphasize the role of matric suction in contributing to the total shear strength of unsaturated soils. The cohesion intercept (c i ) is often defined by:

Critical shear strength equation
Several critical state shear strength equations have also been proposed for unsaturated soils (Alonso et al. 1990;Wheeler and Sivakumar 1995;Wang et al. 2002;Tarantino 2007;Patil et al. 2017). The proposed equations attempt to describe the shear strength of unsaturated soils under critical state conditions in terms of q-p-r space. A general equation form for incorporating matric suction into the critical state shear strength equation can be written as follows: where M is the slope of the critical state line for saturated soils or critical state stress ratio, p is the net mean stress, q is the deviatoric stress. Bishop et al. (1960) proposed a shear strength equation for unsaturated soils as follows:

Discussion of Bishop's parameter
where χ is an empirical parameter. The value of χ has a magnitude between 1 and 0. χ = 1 represents the completely saturated condition, and χ = 0 shows a dry state. However, χ depends on many factors such as soil structure, suction, air-entry value, making it difficult to be estimated. By comparing Eqs. (24) and (34), the value of Bishop's parameter χ can be expressed as below: In this expression, Bishop's parameter χ depends on several factors such as porosity of soils, volumetric water content, residual degree of saturation, disturbance function, which is in good agreement with previous studies (Rojas 2008;Niu et al. 2021). This expression eliminates the necessity of assuming the parameter χ equal to saturation degree solely.

Review of several existing shear strength equations
Numerous equations have been proposed for describing the shear strength of unsaturated soils, which are subdivided into two main categories: fitting equations and estimation equations. Fitting equations are the ones that can be the best fit for a dataset for the determination of one or more fitting parameters, while estimation equations are the ones based on saturated shear strength parameters and additional information. Some equations either directly or indirectly make use of the properties of the SWCC and others utilize additional information such as soil classification or the shear strength at residual suction conditions. It should be noted that most of the equations are empirical and phenomenological. These equations assume that the shear strength of unsaturated soils can be viewed as an extension of saturated shear strength properties. However, the important difference between existing shear strength equations is in the soil property function (SPF), where a different form of SPF will produce a new shear strength equation.
Considering the characteristics of SPF among existing equations, it can be classified into five main equation groups as follows: (1) shear strength equations with SPF depending on volumetric water content (Lamborn 1986;Aubeny and Lytton 2003), (2) shear strength equations with SPF depending on normalized water content (Vanapalli et al. 1996;Fredlund et al. 1996 (5) shear strength equations with SPF depending on residual suction (Rassam and Cook 2002;Naghadeh and Toker 2019). In this study, a representative equation of each group is selected for comparison, and the following section briefly describes the characteristics of the selected shear strength models. Lamborn (1986) proposed a model for predicting the shear strength of unsaturated soils in which the SPF is assumed to be equal to the volumetric water content ( ). This shear strength equation is expressed as follows:

Shear strength equation of Vanapalli et al. (1996)
One of the well-known shear strength equations is the model proposed by Vanapalli et al. (1996), where the SPF is assumed to be equal to the normalized volumetric water content or the so-called effective degree of saturation. The unsaturated shear strength equation is expressed as follows: where θ r = volumetric water content at residual suction, θ s = volumetric water content at saturated condition.

Shear strength equation of Oberg and Sallfors (1997)
Oberg and Sallfors (1997) proposed the SPF to be equal to the saturation degree. The expression for this shear strength model is:

Shear strength equation of Khalili and Khabbaz (1998)
Khallili and Khabbaz (1998)  where (u a − u w ) r is residual suction (or transition suction) that corresponds to residual saturation degree.

Soil-water characteristic curve and evaluation criterion
The proposed model is validated against some of the existing shear strength equations and measured data for six different � soil types: silty clay, sand, silty soil, clay, expansive soil, and sand-kaolin mixture. It should be noted that Eq. (28) of the proposed model is used to obtain results for comparison in the following sections. To predict the unsaturated shear strength of unsaturated soils, the soil-water characteristic curve (SWCC), which shows the relationship between measured suction and degree of saturation, must be presented. In this study, the SWCC model of Fredlund and Xing (1994) is used to plot the reference curve. This model can be expressed as follows: where a, m, n are fitting parameters.
It should be noted that the SWCC and shearing tests are performed on different soil specimens which usually do not have the same boundary conditions (e.g., confining pressure). As a result, the initial void ratios of specimens in those two tests may differ. It is worth noting that a difference in void ratios might produce a change in SWCC, affecting the accuracy of shear strength estimation (Ng and Pang 2000;Lee et al. 2005;Krisdani et al. 2008;Dastjerdi et al. 2014;Oh and Lu 2014;Roy and Rajesh 2018;Zhai et al. 2020;Pham 2022). The calibration procedure for SWCC, therefore, must be undertaken based on the real void ratio of shear test samples in order to overcome the uncertainty of SWCC related to the test conditions. Pham and Sutman (2022) proposed a simplified model for predicting SWCC change with initial density. There are three major advantages to this model: (i) it directly correlates the change in suction with the void ratio, (ii) it predicts the influence of volume change on SWCC using the initial void ratio rather than the current void ratio, and (iii) particle shape and pore size are taken into account within the model. As a result, if SWCC calibration is required due to diverse void ratios, it is recommended to choose this model for implementation, the expression of which is as follows: where e 0 is the initial void ratio of soil sample in SWCC test, which is referred to as reference initial void ratio, e i is the initial void ratio of soil sample in shear test, (u a − u w ) e 0 and (u a − u w ) e i are matric suctions corresponding to the initial void ratios of e 0 and e i , is calibrated constant considering the effect of actual particle shape and size. To determine the parameter , at least two SWCC test data sets at different void ratios are required. When the data of the SWCC test corresponding to different void ratios is unavailable, can be approximated to be between 10 and 15 for sands and silts, between 25 and 35 for clays, and between 45 and 55 for compacted bentonite (Pham and Sutman 2022). Figure 4 shows an example of the SWCC calibration technique utilizing a measured data source from Habasimbi and Nishimura (2019). Figure 4a shows the compression curve where void ratio versus confining pressure (e-lnσ'). It can be observed that the void ratio of the soil sample decreases with increasing confining pressure. For example, the void ratio reduces from 0.726 to 0.711 as the confining pressure is increased from 20 to 600 kPa. Figure 4b shows the SWCC test results at 20 kPa confining pressure (e 0 = 0.726) and calibrated SWCCs at various confining pressures using Eq. (42). It should be emphasized that in Eq. (42), the initial void ratio, e 0 = 0.726, is kept constant as a baseline, and then, the calibrated SWCCs were plotted by replacing the arbitrary void ratio corresponding to different confining pressures. As a result, the effect of volume change and confining pressure on the SWCCs may be adequately modeled. If the confining pressure applied during the SWCC test differs from one of the shear tests, the calibrated SWCCs are recommended to be used to estimate the shear strength of unsaturated soil instead of the tested SWCC results (Pham 2022).
The SWCC results of six different soil types, in which the procedure to determine the air-entry value and residual value are illustrated, are presented in Fig. 5. The physical properties of tested soils and input parameters for shear strength equations are summarized in Table 1.
The performance of shear strength models is assessed through the agreement degree of the curve, which represents the difference between the predicted and the measured curves. The agreement degree of the curve is specialized by average relative error (ARE), which is defined as the percentage of a discrepancy between the value predicted by the analytical equations and the measured one, according to the following expression: where τ measured is the measured shear strength value of ith data, τ predicted is the predicted value of shear strength of ith data, and N is the total number of data points available. Cunningham et al. (2003) conducted a series of triaxial compression tests on the reconstituted silty clay, which was made up of 20% pure Speswhite kaolin, 10% London clay, and 70% silica silt, at varied confining pressures. The slurry soil was pre-consolidated to 130 kPa in an isotropic manner. All of the soil samples utilized in this study were reconstituted soils made from a slurry made at 1.5 times the soil's liquid limit. In a 204 mm diameter lever arm oedometer, the slurry was pre-consolidated one-dimensionally to a maximum vertical effective stress of 200 kPa. After that, the sample was

Comparison results for reconstituted silty clay
discharged under completely drained conditions to ensure that the initial suction was minimal. The interesting point of these test sets is that an air-circulation system was employed to remove moisture from the base of the triaxial sample, and the suctions are independently measured using a samplemounted suction probe. The findings predicted by the proposed model and test data for the suction range of 400 to 1000 kPa are shown in Fig. 6a. When the value ARE is just 1.1%, the proposed model exhibits excellent agreement with measured data for four distinct suction levels. In Fig. 6b, a comparison of shear strength equations and measured data with variations in matric suction is shown, with a constant net normal stress of 400 kPa. The shear strength increases with increasing suction, according to both measured and predicted data. The performance of estimation equations, on the other hand, is vastly different. In comparison to other equations, the proposed model has the best agreement with the measured data. The model of Khalili and Khabbaz (1998) with a value ARE of 26.1% over-predicts significantly while the models of Lamborn (1986) with a value ARE of 50.5% underpredict strongly the unsaturated shear strength. Another finding is that the difference in predicted results between the models of Vanapalli et al. (1996), Oberg and Sallfors (1997), and Naghadeh and Toker (2019) is negligible. This is because the saturation degree of soils, in this case, was relatively high, and the normalized volumetric water contents are approximately equal to the saturation degree. Value ARE of these three models is 13.4%, 12.2%, and 11.3%, respectively. Donald (1956) conducted direct shear tests to measure unsaturated shear strength on well-graded Frankston sand. The Frankston sand specimens were prepared by the vibration technique to obtain the controlled medium dense sand. Once the specimen was prepared on the ceramic disk and set up in the apparatus, it was saturated allowing water, which was in the tank under atmospheric pressure conditions (zero pressure), to flow into the specimen through the ceramic disk. The saturation of the specimen was assumed to have occurred when the change in the mass of the water tank was negligible. At the same time, the normal stress was gradually increased to a predetermined value of 10 kPa. Figure 7 demonstrates the comparison between predicted and measured results for Frankston sands. The suction of sands generally is quite low compared to clays. It is interesting to observe from experimental results that the unsaturated shear strength increases with increasing suction but begins to decrease when suction is continuously increased. Only the proposed model and the one of Vanapalli et al. (1996) is able to reflect reasonably this tendency while other models cannot consider this change. This is because both the proposed model and Vanapalli et al. (1996) model considered the suction contribution to the shear strength by using effective saturation degree. When suction is continuously increased over the residual value, the effective saturation degree becomes very small. As a result, the unsaturated shear strength is reduced although suction increases. However, the proposed model gives a much better agreement with measured data compared to the model of Vanapalli et al. (1996). Value ARE was 1.2%, and 12% corresponding to the proposed model and model of Vanapalli et al. (1996). Additionally, it is found that the model of Khalili and Khabbaz (1998) only matches well with measured data for suction range lower than air-entry value but overpredicts considerably when suction is increased (ARE = 5.7%). It is also observed that the prediction performance of the Lamborn model (1986), and Naghadeh and Toker model (2019) is much lower than other models. The value ARE of these two models is 15.6% and 19.5%, respectively. Gao et al. (2021) conducted a series of suction-controlled triaxial compression tests on compacted silty soils. Silty clay specimens were prepared at the optimum water content (21%) and a dry density of 1.60 g/cm 3 . Firstly, the ovendried kaolin powder and fine quartz sand with a ratio of 3:1 were mixed, and then, distilled water was added to produce the initial water content of 21%. The wet soil was then packed into a cylindrical mold in three layers of roughly equal thickness, and the static compaction method was used to achieve the necessary dry density. Finally, the soil specimen was gently removed from the mold once the compaction process was done. Before the triaxial tests, all specimens were saturated in a vacuum. The axis-translation technique is used to control the matric suction in these data sets. A saturated ceramic disc with a high air-entry value was used to apply pore air pressure to the top of the specimen, while a saturated ceramic disc with a high air-entry value was used to apply pore water pressure to the bottom. A comparison between predicted and measured results for silty soils with the net normal stress of 100 kPa is shown in Fig. 8. It is observed that the results predicted the proposed model is in good agreement with measured data and generally better than other models, with the value ARE of only 1.2%. It is also interesting to note that three models of Vanapalli et al. (1996), Oberg and Sallfors (1997), Naghadeh and Toker (2019) produce a quite close prediction and agree relatively good with measured data. The value ARE for these three models is 2.1%, 2.3%, and 1.8% respectively. The low suction range combined with a high saturation degree can be considered an explanation for the small difference between the three models. It should also be noted that the model of Naghadeh and Toker (2019) is an empirical model that was established based on the test results of silt and the performance of this model seems to be better for silty soils than other soils. On the other side, it is also observed that the model of Lamborn (1986) produces a significant underprediction while the model of Khalili and Khabbaz (1998) gives a high overprediction. Value ARE of these two models are 19.6% and 14.2%, respectively. Kayadelen et al. (2007) investigated the influence of suction on the shear strength of soils using 12 unsaturated triaxial tests. By applying the axis translation technique, the matric suction in the soil specimen may be controlled. The soil samples were saturated prior to the experiments. Using the axis translation, a known magnitude of air pressure was delivered to soil specimens after the saturation process. The soil specimens were then desaturated by eliminating pore water before shearing, and a certain matric suction value was required. Figure 9a presents the comparison between the proposed model and experimental data for a suction range of 0 to 400 kPa. It can be seen that the proposed model produces a good match to measured data for different suction values. The value ARE of the proposed model is only 4.65%. The proposed model is then compared with existing shear strength and measure data for net normal stress case of 50 kPa, and results are indicated in Fig. 9b. It is worthy to note that both analytical and experimental models agree that the shear strength increases nonlinearly with an increase in matric suction. However, the proposed model shows a prediction performance better than other models. The model of Naghadeh and Toker (2019) also produces a good prediction in this case with a value ARE of 7.2% while the value ARE was 11.3% for the Vanapalli et al. model (1996) and 12.9% for the model of Oberg and Sallfors (1997), whereas a low prediction performance is observed for the models of Lamborn (1986), Khalili and Khabbaz (1998), in which the value ARE of these two models is 22.8%, and 20.6%, respectively.

Comparison results for expansive soils
To determine the shear strength of unsaturated expansive soil, Miao et al. (2002) conducted a series of unsaturated triaxial experiments. The samples are remolded expansive soils that have been statically compacted to measure water  Figure 10a shows the comparison between the proposed model and measured data for expansive soils. It is relatively interesting to note that the proposed model agrees well with measured data for all different cases, where variation in matric suction and net normal stress occurs simultaneously. The value ARE of the proposed model, in this case, was only 0.8%. Comparison results between analytical and experimental models are presented in Fig. 10b. It is noted that the model of Khalili and Khabbaz (1998), and the one of Naghadeh and Toker (2019) produce a high overprediction of unsaturated shear strength. The value ARE of the two models is 35% and 16%, respectively. On the other hand, the model of Lamborn (1986) highly underpredicts the shear strength of unsaturated soils, which is specialized by the value ARE of 26.2%. It is also observed that the prediction performance of models of Vanapalli et al. (1996), andOberg andSallfors (1997) decreases with increasing suction. The value ARE was 7.53% for the model of Vanapalli et al. (1996), and 8.33% for Oberg and Sallfors (1997). Guan et al. (2010) conducted a series of shear strength tests on the sand-kaolin mixture using modified triaxial apparatus.  Vanapalli et al. (1996) Oberg & Sallfors (1997) Khalili & Khabbaz (1998) Naghadeh & Toker (2019) The sand-kaolin mixture comprises 35% of sand and 65% of kaolin. The sand-kaolin specimens were statically compacted to obtain the maximum dry density of 1.67 g/cm 3 and the optimum water content of 19%. All specimens were saturated at the beginning of the test to have a uniform initial condition and to reduce the matric suction to a low value. Saturation was performed by applying cell pressure and back pressure from digital pressure and volume controller equipment. The specimen was then isotropically consolidated to the designated net confining pressure. Once consolidation was completed, the drying process of the specimen was conducted by applying the designated matric suction using the axis translation technique. Figure 11 shows a comparison between predicted and measured results for the sand-kaolin mixture. Like other soil types, it can be observed that the model of Khalili and Khabbaz (1998) produces a high overprediction while the model of Lamborn (1986) provides a significant underprediction of unsaturated shear strength. The value ARE of the two models was 19.3% and 25%, respectively. The models of Vanapalli et al. (1996), Oberg andSallfors (1997), and Naghadeh and Toker (2019) also underpredict slightly the shear strength of unsaturated sand-kaolin mixture. The value ARE of the three models was 9.02%, 9.76%, and 5.82%, respectively. However, it should also be noted that the performance of these three estimation models tends to decrease with an increase in matric suction. It is observed that the proposed model gives a good agreement with measured data and is generally better than other models when the value ARE was only 2.1%.

Statistical evaluation
A statistical evaluation was conducted on six different soil types for shear strength models, as shown in Fig. 12. It is interesting to note that the value ARE of the proposed model was less than 10% among all six cases which reveals that the proposed equation has a good and reliable performance in predicting the shear strength of unsaturated soils. It is also noted that the models of Lamborn (1986), Khalili and Khabbaz (1998) show a lower performance for different soil types, among the considered models. Definition of SPF based on volumetric water content or air-entry value can be considered the reason behind the lower performance. It should be noted that the prediction performance of models of Vanapalli (1996), Oberg and Sallfors (1997), and Naghadeh and Toker (2019) is quite similar. However, the model of Vanapalli (1996), and the model of Naghadeh and Toker (2019) are much sensitive to the accuracy of residual value. Another interesting finding is that the prediction performance of shear strength models for the clayey soils and expansive soils is lower than that for silty and sandy soils.
The measured against predicted shear strength is presented in Fig. 13. The results indicate that the proposed model shows an excellent performance in estimating shear strength for any suction range. It is also found that the existing shear strength models tend to agree well with measured data for low shear strength range but differ significantly  Vanapalli et al. (1996) Oberg & Sallfors (1997) Khalili & Khabbaz (1998) Naghadeh & Toker (2019) when the shear strength range is increased. It is interesting to note that the model defining SPF based on residual suction value (Naghadeh and Toker 2019) gives a better prediction than the one defining SPF based on air-entry value (Khalili and Khabbaz 1998).

Non-isothermal shear strength model for unsaturated soils
It is well known that climatic conditions subject soils near the ground surface to cyclic heating and cooling episodes. However, the shear strength characteristics of soils are highly affected by temperature variations. As temperature increases, the reduction of matric suction and void ratio triggers the change of shear strength. Therefore, it is a necessity to estimate the shear strength of soils with temperature change. As a result, the proposed model is further extended to predict the effect of temperature on the unsaturated shear strength of soils in this section. It should be noted that a temperature increase tends to induce a variation in the shear strength of soils through four components: friction angle, cohesion, disturbance function, and matric suction. For the analysis, the following notation will be used: a subscript T indicates the current temperature, and a subscript T 0 indicates the reference temperature. The Matric suction, u a -u w (kPa) Measured data Proposed model Lamborn (1986) Vanapalli et al. (1996 Oberg & Sallfors (1997) Khalili & Khabbaz (1998) Naghadeh & Toker (2019) reference temperature is defined as the initial temperature of the soil sample during matric suction measurement in the laboratory and is often the same as room temperature. Based on Eq. (28), the non-isothermal shear strength model for unsaturated soils can be expressed as follows: where T = unsaturated shear strength at temperature T, D T is disturbance function at temperature T, u a − u w T is matric suction at temperature T, c ′ T and ′ T is cohesion and friction angle at temperature T. However, many researchers (44) found that the parameters c ′ T and ′ T are either independent of temperature change or the impact of temperature change is limited in most cases (Hueckel and Baldi 1990;Graham et al. 2001;Cekerevac and Laloui 2004;Laloui and Sutman 2021).

Temperature-dependent matric suction
In porous media, the mechanical and thermodynamic equilibrium is often derived from the relationship between capillary potential and the free energy of soil moisture (Arya and Paris 1981;Pham and Suman 2022). The matric suction or capillary pressure is determined by the following relation:  where s is surface tension, cosα is contact angle, r i is pore radius, R is particle radius, e is the void ratio, n s is the number of particles for a soil volume, δ is a calibrated coefficient considering the real shape and orientation of soil particles.
The change of matric suction with temperature thus can be specialized through variation in surface tension, contact angle, particle radius, and void ratio as follows: The matric suction at the current temperature u a − u w T can be related to the reference matric suction u a − u w T 0 by the simplified form as follows: in which f is the temperature-dependent function of surface tension, f is the temperature-dependent function of air-water contact angle, f R is the temperature-dependent function of particle radius, f e is the temperature-dependent function of void ratio, v is volumetric thermal expansion coefficient of solids, e 0 and e T are initial void ratios at reference and current temperature respectively. It should be noted that T is expressed in degrees Kelvin (K) in the above formulae.

Temperature-dependent disturbance function
The disturbance function is known to depend on the void ratio of soils. The void ratio of soils generally varies with the temperature change. The disturbance function considering thermal volume change, therefore, can be re-written as follows: where D 0 is the disturbance function at reference temperature (T 0 ). Demars and Charles (1982) found that the void ratio variation of normally consolidated soils due to temperature  Vanapalli et al. (1996) Oberg & Sallfors (1997) Khalili & Khabbaz (1998) Naghadeh & Toker (2019) Linear line 1: 1 fluctuation is directly related to soil plasticity. Based on experimental outcomes, an equation that describes the relationship between change in void ratio due to temperature cycle and plasticity index is proposed as follows: where I P is the index of plasticity, ΔT is temperature increment (°K).

Validation of non-isothermal shear strength model
The validation of the non-isothermal shear strength model is investigated against experimental results on the unsaturated compacted silt. Uchaipichat and Khalili (2009) (53) e T = e 0 − 0.00048 + 0.0000088I P ΔT conducted a comprehensive program of non-isothermal testing on a compacted sample of silt in a modified triaxial cell. Applied temperature values ranged from 25 to 60 °C, and the suction values varied from 0 to 300 kPa. The soil was statically compacted to a dry unit weight of about 1.53 g/cm 3 , and at a moisture content of 10.5%. The samples were prepared dry of optimum to render the soil matrix amenable to stiffening with increasing matric suction. Before compaction, the soil was air-dried at room temperature, ground by a rubber hammer, and screened through a 400 μm sieve. Then, it was carefully wetted with a spray gun to a water content of 10.5% and placed in a sealed plastic bag to cure for 24 h for moisture equalization. It should be noted that the variation in soil-water characteristic curve with temperature (non-isothermal SWCC) can be predicted by combining Eqs. (41) and (47) in this study. Figure 14 shows the comparison between predicted and measured results for the non-isothermal SWCC. It is observed that the predicted SWCC matches well with measured data for different values of temperature. Both analytical and experimental models agree that higher temperature produces lower SWCC, which corresponds to a reduction in suction with increasing temperature. However, the effect of temperature on SWCC is only significant for the suction range in the transition zone (between air-entry value and residual value). Figure 15 presents comparison outcomes for unsaturated shear strength at elevated temperatures. It is interesting to note that the unsaturated shear strength decreases with increasing temperature. However, the influence degree of temperature on shear strength depends strongly on the suction range. For example, the unsaturated shear strength decreases 14.2% at the suction of 100 kPa and 18.6% at the suction of 300 kPa when the temperature is increased from 25 to 60 °C. It is also worthy to note that the value ARE was 4.2%, which indicates that the proposed model produces a high satisfaction with measured data for predicting the variation in shear strength with temperature.

Conclusions
In this paper, the evolution of the shear strength of unsaturated soils concerning the matric suction was studied under isothermal and non-isothermal conditions. Several main conclusions can be drawn from this study as follows: A new shear strength equation based on the micromechanical model and the disturbed state concept was presented for unsaturated soils. The original point of this study is considering the solid contact area ratio which was neglected among most current models. Using the proposed model, the non-linear relationship between the matric suction, saturation degree, and the shear strength of unsaturated soils is described.
The results obtained from the proposed model were compared to several existing shear strength models and measured data for six different types of soils. The results indicate that the proposed model has a good performance in predicting the shear strength of unsaturated soils, and generally is better than other existing models. It is also found that the models defining SPF based on effective saturation degree give a better prediction than the models defining SPF based on the air-entry value or volumetric water content.
In response to varying climatic conditions, the proposed model was then extended to produce a non-isothermal shear strength equation for unsaturated soils. Validation of non-isothermal shear strength and SWCC models was performed against the experimental data on compacted silt. The comparison results show a good performance of the proposed model in predicting the shear strength variation with temperature.
Notation A: Total section area (m 2 ); A c : Solid contact area between the two soil grains (m 2 ); A w : Total water area (m 2 ); e: Void ratio (dimensionless); f: Water distribution function (dimensionless); f : Temperature-dependent function of surface tension (dimensionless); f : Temperature-dependent function of air-water contact angle (dimensionless); f r : Temperature-dependent function of particle radius (dimensionless); f e : Temperaturedependent function of void ratio (dimensionless); : Total stress (Pa); c : Solid contact stress (Pa); ′ : Effective stress (Pa); c ′ : Effective cohesion (Pa); c ′ T : Cohesion at temperature T (Pa); c i : Cohesion intercept (Pa); D T : Disturbance function at temperature T; D: Disturbance function or solid contact area ratio (dimensionless); I P : Index of plasticity (dimensionless); M: Slope of the critical state line (dimensionless); N: Total number of data points available; n: Porosity of soils (dimensionless); n s : Number of contact points (dimensionless); P: Total external load applied on the unsaturated soil (N); P c : Force transferred through the grain contact (N); P w : Total force carried by pore water (N); P a : Force carried by pore air (N); P m w : Forces carried by water in two menisci (N); p: Net mean stress (Pa); q: Deviatoric stress (Pa); u w : Pore water pressure (Pa); u a : Pore air pressure (Pa); u a − u w : Matric suction (Pa); (u a − u w ) r : Residual suction (Pa); u a − u w T : Matric suction at temperature T (Pa); R : Particle radius (m); S e : Effective saturation degree (dimensionless); S r : Residual degree of saturation (dimensionless); T: Current temperature (Kelvin degree); t: Contact thickness (m); V f w : Free water volume (m 3 ); V m w : Menisci water volume (m 3 ); V c : Total solid contact volume (m 3 ); V ci : Contact volume between two solids (m 3 ); V s : Solid volume (m 3 ); V si : Volume of a solid particle (m 3 ); ′ : Friction angle of soils at saturated condition (degree); ′ T : Friction angle at temperature T (degree); τ measured : Measured shear strength value of ith data (Pa); τ predicted : Predicted value of shear strength of ith data (Pa); T : Unsaturated shear strength at temperature T (Pa); : Contact angle (radian); χ: Bishop's parameter (Pa); e : Normalized volumetric water content (dimensionless); : Volumetric water content (dimensionless); θ r : Volumetric water content at residual suction (dimensionless); θ s : Volumetric water content at saturated condition (dimensionless); v : Volumetric thermal expansion coefficient of solids (dimensionless) Abbreviations DEM: Discrete element method; SWCC : Soil-water characteristic curve; AEV: Air-entry value; ARE: Average relative error; SPF: Soil property function Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.