Bounding surface constitutive model for unsaturated soils considering microscopic pore structure and bonding effect

. This article presents a bounding surface model for unsaturated soils by using skeleton stress and bonding variable based on microcosmic pore structure as constitutive variables. A Hydraulic hysteresis soil-water characteristic curve model considering deformation and hydraulic hysteresis is combined to achieve hydraulic coupling. The proposed model can capture the change of the inter-particles bonding effect in the deformation process of unsaturated soils and accurately predict the hydraulic mechanical behavior of unsaturated soils under complicated loading paths and wetting-drying cycles. The validity of the proposed model is confirmed by the results of unsaturated isotropic compression tests, wetting-drying cycles tests and unsaturated triaxial shear tests reported in the literature.


Introduction
Unsaturated soil has a bonding effect due to the existence of water meniscus between particles, which makes the mechanical properties of unsaturated soils more complicated than saturated soils. It is of great value to establish unsaturated constitutive models considering the bonding effect. The bonding effect is related to suction, degree of saturation, and pore structure. Some researchers indirectly describe the bonding effect by considering the influence of suction and saturation in their model [1][2][3][4][5][6]. Gallipoli et al [7] were the first to establish an elastic-plastic constitutive model that directly describes the bonding effect of unsaturated soils by using a bonding variable to replace the suction.
Besides suction and degree of saturation, the effect of microscopic pore structure on the hydraulic properties of unsaturated soils are not negligible. Some researchers [8][9][10] found that the pore structure of unsaturated soils has a significant influence on the bonding effect and hydraulic mechanical behavior between particles. But the existing models have little consideration for the effect of microscopic pore structure. In addtiton, compared with the traditional elasto-plastic model, the bounding surface model [11] is versatile and can more accurately predict the mechanical behavior of various soil types [12][13]. For unsaturated soils with more complicated mechanical properties, it is more suitable to use the bounding surface plastic model [14][15].
In this paper, a hydraulic coupling bounding surface model for unsaturated soils using improved average skeleton stress and bonding variable as basic constitutive variables is established. Experimental results of isotropic compression tests, wetting-drying cycles tests and triaxial shear tests are used to validity the model.

Constitutive Variables of Model
The improved average skeleton stress and bonding variable are used as the constitutive variables of the model in this paper to consider the influence of microscopic pore structure and bonding effect. Based on the research results of Alonso et al [16], the parameter χ of Bishop effective stress [17] is replaced by the effective degree of saturation S e . The improved average skeleton stress is written as follows: Where p net is the net stress, S e is the effective degree of saturation that can be express as Where S r is the degree of saturation, S res is the residual degree of saturation The degree of saturation S r of the bonding variable defined by Gallipoli et al is replaced by the effective degree of saturation S e to establish a new bonding variable as follows : Where the f(s) is the ratio of the stabilising interparticle force at the two suctions of s and zero for the ideal case of a water meniscus located at the contact between two identical spheres; (1−S e ) is the accounts for the number of water menisci between the macropores per unit volume of solid fraction.

Unsaturated Normally Consolidated Lines
Establishing the relationship between unsaturated soil deformation and bonding variable is important to determine the normally consolidated lines for unsaturated soils. Using a methodology similar to Gallipoli et al [7], the authors assume that a unique functional relationship is satisfied between e/e s and ζ : Where e and e s are the void ratio corresponding to the unsaturated and saturated at the same value of improved average skeleton stress p * , a and b are fitting parameters, e s can be calculated according to the saturated normal compression line: Where N is the void ratio corresponding to p * = 1; λ is the slope of the saturated normal compression line in the ln p * − e s plane.
Thus, the unsaturated normal compression curve in the ln p * − e s plane can be calculated by combining Eq.

Yield equation
Eq. (6) defines a normal compression surface in the ln p * − e −ζ space. The irreversible volumetric strains occur when the stress path on the normal compression. Yield equation can be express as follows according to Hu et al [18] ( )

Bounding surface
The yield surface of Modified Cam-clay model [19] is used to be the bounding surface in the proposed model. Fig. 1 shows the bounding surface in three-dimensional state, the bounding surface can be expressed as follows: Wherep * is the improved average skeleton stress on the bounding surface,q is the deviator stress on the bounding surface,p c * [ζ, p c * (0)] is the yield stress corresponding to constant value of bonding variable ζ defined by Eq. (7) , M is the slope of the critical state line.

Flow rule
A non-associated flow rule is chosen in the proposed model. The plastic potential surface equation is written as follows: Where α is a parameter determined by zero lateral strain when the stress state corresponds to the K 0 value [20], according to the derivation of Alonso et al [1], α can be expressed as The incremental forms of plastic strain and plastic shear strain can be obtained as follows: Where Λ is the plastic multiplier. According to the plastic theory of bounding surface, plastic multiplier Λ can be expressed as follows: WhereK p is the plastic modulus corresponding to the mapping point on the bounding surface.

Hardening law
The yield stress p c * (0) is chosen as the hardening parameter of the proposed model. The hardening law can be obtained as follows: Where dε v p is the plastic volumetric strain increment.
The plastic modulus on the bounding surface can be obtained by combining the flow rule, hardening law and the condition of consistency on the bounding surface.

Mapping rule
Zienkiewicz et al [21] establishes radial mapping rules with simple form, few parameters, and good applicability. A similar radial mapping rule is used in the proposed model: where δ is the distance between the current stress point and the mapping origin, δ 0 is the distance between the mapping point and the mapping origin, r is the mapping exponent related to the basic physical properties of the soils: where r 0 is an initial value, n is proportionally constant (n＞0), p dε  is cumulative value of plastic strain increment.
As shown in Fig. 2, the mapping origin is the coordinate origin in the p * − q plane, and the mapping point on the bounding surface is determined by the intersection of the line connecting the mapping origin to the current stress point and the bounding surface.
where G is elastic shear modulus.

Hydraulic behaviour
Gallipoli et al [22] developed a soil-water characteristic curve model considering deformation effect based on the van Genuchten model [23], which is very representive: Where v is specific volume, m 1 , m 2 , m 3 m 4 are fitting parameters.
Based on Gallipoli et al [22], the authors establish a soil-water characteristic curve model considering hydraulic hysteresis to achieve hydraulic coupling. The main drying curve and the main wetting curve are defined as follows: Where m 1d is fitting parameters for the main drying curve, m 1w is fitting parameters for the main wetting curve.
The scanning curve equations are expressed as follow based on Zhou et al [6] ( ) where S ed is the effective degree of saturation corresponding to the main drying curve, S ew is the effective degree of saturation corresponding to the main wetting curve, s d is the suction corresponding to the main drying curve at the same effective degree of saturation as the current point in the s−S e plane, s w is the suction corresponding to the main wetting curve at the same effective degree of saturation as the current point in the s−S e plane, η is the parameter that controls the shape of the scanning curves.
Where S eds is the effective degree of saturation corresponding to the drying scanning curve, S ews is the effective degree of saturation corresponding to the wetting scanning curve.

Model verification
A series of experimental data of reconstituted kaolin [24], mixture of bentonite and kaolin [25] Fig. 3, Fig. 4, Fig. 5 show the comparisons between model predictions and experimental data.

Isotropic compression tests at constant suction
Sivakumar [24] conducted a series of isotropic compression tests of reconstituted kaolin at constant suctions. Two sets of tests at suctions of 200 kPa, and 300 kPa are used to verify the proposed model. The isotropic loading paths experienced by the samples under each suction is 50 kPa→200 kPa, 50 kPa→250 kPa. Fig. 4 shows the comparisons between model predictions and experimental data. The model can accurately predict the variation of void ratio. As the net stress increases, the void ratio gradually decreases.

Wetting-drying cycle test at constant net stress
The experiment results of constant net stress wettingdrying cycle tests conducted by Sharma [25] are used to verify the performance of the proposed model. The suction path is 300 kPa→20 kPa→300 kPa. Fig. 5 shows the comparisons between experimental results and model predictions during wetting-drying cycle. It can be seen that the proposed model can effectively capture the hydraulic hysteresis characteristics of unsaturated soils during wetting-drying cycle.

Triaxial shear tests at constant suction
Sivakumar [25] conducted a series of triaxial shear tests on reconstituted kaolin at the suctions of 200 kPa, 300 kPa. The samples first experience a consolidation with a net stress path of 50 kPa→100 kPa and then maintain the stress path of Δq/Δp net =3 for drainage shearing.
The comparisons between experimental results and model predictions of constant suction triaxial shear tests are shown in Fig. 6. Clearly the proposed model can effectively predict the mechanical behavior of unsaturated soil experiencing triaxial shearing. The deviator stress increase with the increase of the axial strain, the specific volume decrease with the increase of the axial strain, finally reaches the critical state. As the suction increases, the deviator stress increases gradually.

Conclusions
A new hydraulic coupling bounding surface model for unsaturated soils considering microscopic pore structure and bonding effect is established in this article. Using the effective degree of saturation to replace the degree of saturation, a new average skeleton stress and bonding variable are established as the constitutive variables of the model, the plastic deformation is calculated by the bounding surface plasticity theory, and the hydraulic coupling is achieved by combining the hydraulic hysteresis soil-water characteristic curve model considering the influence of deformation. The proposed model can effectively predict the mechanical behavior of unsaturated soils under isotropic loading, shear loading and wetting-drying paths.