Hydraulic hysteresis effects on the coupled flow–deformation processes in unsaturated soils: Numerical formulation and slope stability analysis

Accepted Manuscript

Hydraulic hysteresis effects on the coupled flow-deformation processes in unsaturated soils: Numerical formulation and slope stability analysis Ran Hu , Jia-Min Hong , Yi-Feng Chen , Chuang-Bing Zhou PII: DOI: Reference:

S0307-904X(17)30574-7 10.1016/j.apm.2017.09.023 APM 11968

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

30 March 2017 11 August 2017 6 September 2017

Please cite this article as: Ran Hu , Jia-Min Hong , Yi-Feng Chen , Chuang-Bing Zhou , Hydraulic hysteresis effects on the coupled flow-deformation processes in unsaturated soils: Numerical formulation and slope stability analysis, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.09.023

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights A coupled flow-deformation numerical method with hydraulic hysteresis is proposed.



A return mapping algorithm is developed to implement the hydraulic hysteresis.



The performance and convergence of the numerical formulation is improved.



The significant effects of hydraulic hysteresis on slope stability is demonstrated.

AC

CE

PT

ED

M

AN US

CR IP T



1 / 46

ACCEPTED MANUSCRIPT

Hydraulic hysteresis effects on the coupled flow-deformation processes in unsaturated soils: Numerical formulation and slope stability analysis

a

CR IP T

Ran Hua, Jia-Min Honga,b, Yi-Feng Chena* and Chuang-Bing Zhoua,c

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China b

Powerchina Huadong Engineering Corporation Limited, Hangzhou 310004,China

School of Civil Engineering and Architecture, Nanchang University, Nanchang 330031, China

AN US

c

*

M

Corresponding author: [email protected], [email protected]

Aug 11, 2017

AC

CE

PT

ED

Revised Manuscript

A paper submitted for possible publication in Applied Mathematical Modelling

2 / 46

ACCEPTED MANUSCRIPT

Abstract: The hysteresis of water retention curve has a profound influence on the coupled hydro-mechanical behaviors in unsaturated soils, but numerical implementation with consideration of this property was rarely reported due to the difficulties in the integration of the coupled constitutive models. In this study, a numerical formulation is proposed for modeling the coupled flow-deformation processes with hydraulic hysteresis. A return mapping scheme is

CR IP T

developed to integrate the water retention curve model with hydraulic hysteresis and the elasto-plastic model simultaneously within a time step, and the deformation-dependent nature of the water retention curve is considered rigorously by modifying the coefficient matrices in the discretized governing equations. The performance and efficiency of the proposed numerical formulation is validated by two existing laboratory tests and a computational example,

AN US

demonstrating better performance and convergence of the proposed formulation. The proposed procedure is then applied for modeling the coupled flow-deformation processes in a soil slope under rain infiltration. The simulated results reveal the significant effects of hydraulic hysteresis on the coupled water-air two-phase flow and elasto-plastic deformation processes. The solid

M

deformation and the evolution of the shear band would be remarkably overestimated, and the

PT

ED

slope failure would be early predicted when neglecting hydraulic hysteresis.

CE

Keywords: unsaturated soils; coupled hydro-mechanical analysis; hydraulic hysteresis;

AC

numerical formulation; landslides

3 / 46

ACCEPTED MANUSCRIPT

1 Introduction Understanding and modeling the coupled fluid flow and deformation processes in unsaturated soils is crucial to many engineering applications, including rainfall-induced landslides [1-6], nuclear water disposal [7-10] and embankment/rockfill dam engineering [11-13]. Although it has long been studied, characterizing the coupled flow-deformation processes is still a challenging

CR IP T

issue, due mainly to the strong solid-liquid interactions that greatly complicate the development of constitutive models and the difficulties in implementing the constitutive models into a computer code for the coupled hydro-mechanical analysis. A large number of researchers have focused on the development of constitutive models for reasonably and comprehensively

AN US

describing the coupled stress-strain and hydraulic behaviors of unsaturated soils [14-28]. But much less emphasis has been given to the numerical implementation and application of the coupled hydro-mechanical models [29-42].

One of the challenges involved in the numerical implementation is the integration

M

procedure for the coupled constitutive models with hydraulic hysteresis of the water retention curve. Although the classic implicit and explicit schemes have been widely used [43-45], these

ED

classic integration methods need to be modified for the integration of the coupled hydro-mechanical constitutive models for unsaturated soils, because the stress state and

PT

hardening values for the mechanical behavior as well as the degree of saturation for the hydraulic behavior should be updated simultaneously. More importantly, due to hydraulic

CE

hysteresis and deformation-dependence of water retention behavior, the degree of water

AC

saturation is related to not only the current hydraulic state (denoted by the suction and the void ratio), but also the hydraulic path. The integration of the water retention curve model with hydraulic hysteresis significantly complicates the numerical implementation, and the contributions to this topic can be found in [23, 29, 30, 38, 46-49]. Vaunat et al. [38]

presented

an implicit integration scheme implemented based on the return mapping method [43] for the Barcelona Basic Model [14]. Sheng et al. [30] later presented an explicit stress integration scheme for unsaturated soils with automatic sub-increment and error control method by Sloan [44]. Zhou and Sheng [23] developed an adaptive sub-stepping explicit integration scheme for a 4 / 46

ACCEPTED MANUSCRIPT

porosity-dependent hydro-mechanical model for unsaturated soils. Hu et al. [29] developed a return-mapping implicit integration method for coupled hydro-mechanical constitutive model with a focus on the deformation-dependence of water retention curve. For the above integration schemes, although hydraulic hysteresis (including its dependency on void ratio) was fully considered in a couple of constitutive models [22, 23, 26], the hydraulic path was always assumed

CR IP T

to remain on the main wetting or drying surfaces [5, 12, 29, 50] and very few works focused on the influence of hydraulic hysteresis on the coupled hydromechanical processes [40, 42]. On the other hand, owing to the deformation-dependence of the water retention curve model, the derivation of the governing equations for the two-phase fluid flow and deformation processes can induce an additional term of the first-order partial derivative of saturation with respect to

AN US

void ratio. This term was also neglected in the previous studies [3, 12, 51] and should be rigorously considered.

It has been well understood that hydraulic hysteresis is fundamental to quantifying the soil-water-air coupling in the application of rainfall-induced landslides [4, 5, 50, 52-57]. As the

M

wetting front propagates during rain infiltration, the hydraulic path of soil in slope can remain on

ED

the wetting surface or extend along the scanning curve, depending on the current hydraulic state and its history [32]. Due to the reduction of soil shear strength with strain rate during the wetting

PT

front propagation, single or multiple slip surfaces develop progressively, which refers to progressive failure. In the analysis of progressive failure of slope, hydraulic hysteresis has

CE

always been neglected and the main drying curve has always been employed because the measurement of the main wetting curve is more time-consuming. Some studies have been

AC

conducted to examine the effect of the main wetting/drying surfaces on the slope deformation. Ebel et al. [57] investigated the hysteretic effects on the slope stability and concluded that the potential for landslide in unsaturated materials may be underestimated, which was also reported by Tsai [52]. The experimental observation and numerical simulation by Ma et al. [53] showed that the water retention curve hysteresis substantially affects the water content distribution in, and the stability of, the slope. However, most of these studies focused on the effect of hydraulic path on the main wetting/drying surface, and very few studies reported the influence of hydraulic

5 / 46

ACCEPTED MANUSCRIPT

hysteresis (including main wetting/drying as well as scanning curves) on the slope deformation and stability [32]. This study develops a modeling approach for the coupled flow-deformation analysis with consideration of hydraulic hysteresis of water retention curve, to examine the impact of hydraulic hysteresis on the coupled flow-deformation processes in unsaturated soils. After

CR IP T

briefly reviewing a coupled hydro-mechanical model for unsaturated soil in Section 2, in Section 3, the coupled flow-deformation coefficient matrix representing the deformation-dependent water retention curve is developed, and an implicit integration procedure is then proposed to update the degree of saturation including the main wetting/drying surface and scanning curves. In Section 4, the performance of the proposed numerical model is validated by two laboratory

AN US

tests and a computational case of rain infiltration into a soil column. Finally, the proposed numerical method is applied to quantify the effect of hydraulic hysteresis on the coupled

M

unsaturated flow and elasto-plastic deformation processes in a slope under rain infiltration.

2 Review of a coupled hydro-mechanical constitutive model with

ED

hydraulic hysteresis effects

PT

A coupled hydro-mechanical constitutive model is essential to simultaneously predict the stress-strain and the water retention behaviors in unsaturated soils. Based on the developed

CE

hysteretic water retention curve [58] and the elasto-plastic model [59], a coupled model was formulated for unsaturated soils within the framework of thermodynamic [22], and this coupled

AC

model will be briefly introduced below. By the nature of hysteresis, the degree of saturation is not only related to the current

hydraulic state, but also dependent on the hydraulic path. The water retention behavior can be categorized into two groups, including the main wetting/drying surfaces and the scanning curves. In our previous work, a deformation-dependent hysteretic water retention curve model was developed based on the change in pore size distribution of soils subjected to hydraulic and mechanical loading. As shown in Fig. 1, the water retention curve model defines a bounding 6 / 46

ACCEPTED MANUSCRIPT

surface on which the hydraulic response is either irreversible (denoted by domain DM) or reversible (denoted by domain DS):

DM  ( Se , s, e) | Se =Se, w or Se =Se, d 

(1)

DS  ( Se , s, e) | Se, w  Se  Se, d 

(2)

where Se is the effective degree of saturation, transformed into Sr by Sr  Sr, res + (1  Sr, res)Se, with

CR IP T

Sr, res being the residual degree of saturation.

The main drying (Se, d) and main wetting (Se, w) surfaces given in Eqs. (1) and (2) can be defined by

Se,  [1  (  exp(kpe)s)nw ] mw ,  w, d

(3)

AN US

where s is the suction, e is the void ratio, nw, mw, kp and are model parameters  d denotes the main drying path, and   d denotes the main wetting path.

Any infinitesimal movement of a soil state along the bounding surface DM will induce an irreversible change in saturation, given by

M

ds  dS r  (1  S r,res )Se (1  Se1 mw )  1 s kp    de 

(4)

ED

The following relation is adopted to describe the scanning curves in the reversible domain DS:

PT

ds  dS r  (1  S r,res )Se (1  Se1 mw )   kss s kse    de 

(5)

CE

where kss and kse are, respectively, the slopes of the asymptotes for the main hydraulic and mechanical scanning curves in the lnSelns and lnSee planes (Fig. 1).

AC

It is worth noting that the suction s and the void ratio e in Eqs. (4) and (5) are not independent.

The classic water retention curve was always determined by fitting the suction-saturation data under constant stress conditions during which the void ratio could change due to the evolution of suction, suggesting that the influence of void ratio (or deformation) on the water retention properties should be embedded in the parameters of the classic water retention model. In this study, however, the influence of void ratio (or deformation) on water retention curve is explicitly taken into account by Eqs. (4) and (5). Although it is difficult to perform experiments to directly

7 / 46

ACCEPTED MANUSCRIPT

measure the suction-saturation data under constant void ratio, the parameters of the water retention curve model in this work could be calibrated by the contour map of void ratio in the suction-saturation plane after the data of suction-saturation-void ratio were measured [60, 61]. On the other hand, given the fact that the inter-particle bonding plays a key role in the deformation of unsaturated soils [16, 17], the elasto-plastic model introduces a dimensionless

1  S r1/4  0.32e2  4.06e  0.11 The average skeleton stress is defined as [e.g., 62]:

(6)

(7)

AN US

     pa   Sr s

CR IP T

variable (bonding factor, ) to represent this effect. The bonding factor  was defined as [59]:

where  is the average skeleton stress tensor,  is the total stress tensor,  is Kronecker delta and

pa is the net stress tensor. Eq. (7) has the same expression of Bishop’s stress [63] if the weighting factor, , is replaced by the degree of saturation. Note that the effective stress parameter

M

 can be taken as the degree of saturation only in the condition that the work of water-air interface is neglected [64, 65]. Khalili and Khabbaz [66] and Khalili et al. [67] proposed a relation for this

ED

parameter in terms of the suction ratio s/se, where se is the value of suction corresponding to the transition between saturated and unsaturated states. In this study, although the work of water-air

PT

interface is neglected by employing the average skeleton stress (Eq. (8)), the role of the interface can be taken into account by introducing the bonding factor (Eq. (6)) [22].

AC

CE

By Eq. (6), the loading-collapse curve can be formulated by [59]:

 (   )ln pc (0)  Na b  pc ( )  exp        a b  

(9)

where pc() is the related yield average skeleton stress under isotropic stress conditions for a specific value of , N is the void ratio corresponding to p1 kPa,  is the slope of the saturated normal compression line in the lnp-e plane,  is the slope of the elastic swelling line in the lnp-e plane, pc(0) is the yield stress at saturated state, and a and b are fitting parameters. The yield function (f) and the plastic potential (g) can be developed based on Modified Cam-Clay model [68]: 8 / 46

ACCEPTED MANUSCRIPT

where

f  q2  M 2 p[ pc ( )  p]  0

(10)

g  q2  M 2 p[ pc ( )  p]  0

(11)

p   13  ii is the mean average skeleton stress, q is the deviatoric stress, q  3J 2 , J2is

the second deviatoric stress invariant, M is the slope of the critical state line in qp plane,  is a determined

by

imposing

zero

lateral

strains

M(M9)(M3)/[9(6M)()] [69].

under

K0

conditions,

CR IP T

constant

The incremental plastic strain ( d v ) and plastic deviatoric strain ( d q ) can be obtained by p

p

the flow rule:

AN US

p  g p d v    p    d q    g q  

(12)

where  is a positive plastic multiplier (0).

Fig. 2 presents a three-dimensional view of the initial yield surface and the subsequent yield surface. The hardening rule of the yield surface is expressed as

M

dpc (0) 1  e  d vp pc (0)   

(13)

ED

p where d v is the incremental plastic volumetric strain determined by Eq. (12).

e  0  dp d v   1 K  e     d q    0 1 (3G )   dq  

(14)

CE

PT

The elastic response of the stress-strain behavior is given by

where d v and d qe are the incremental elastic volumetric strain and incremental deviatoric e

AC

strains, respectively, K(1e)p/ and G3K(12)/[2(1)] are the elastic bulk and shear moduli, respectively, and  is Poisson’s ratio. The permeability of the deformable soils can be expressed as [58] k ( e) 

k0  2 exp(2kp e)  exp(2kp e0 ) 2 0

(15)

where e0 is the reference (initial) void ratio, kp is the parameter involved in the water retention curve model given in Eq. (3). 9 / 46

ACCEPTED MANUSCRIPT

The relative permeability functions for the water flow (krw) and air flow (kra) have the following forms [58]:

krw ( Se )  Se1 2 [ I S1/mw ( mw  1 nw ,1  1 nw )]2

(16)

kra ( Se )  (1  Se )1 2 [1  I S1/mw ( mw  1 nw ,1  1 nw )]2

(17)

e

e

z

 t  1  t  I ( x, y )   t  1  t  0 1

z

x 1

y 1

dt

x 1

y 1

dt

0

(18)

AN US

3 Finite element formulation

CR IP T

where Iz(x,y) is the incomplete Beta function,

3.1 Governing equations

Unsaturated soils can be regarded as a three-phase deformable porous medium, i.e. soil skeleton

M

for solid phase (s), water for liquid phase (w) and air for gas phase (a). In this work, we neglect the transfer of mass among the three phases and only consider the small deformation under isothermal

ED

condition. From the mass conservation for the water phase and gas phase, the governing equations for water/gas flow can be written as [62]:

PT

  Sr  B   2     S r  pw  B S r (1  S r ) pa  Ks Ks  Kw 

CE

   kk       B Sr s  S r      rw  pw   w g     B Sr v  0 Ks    w 

AC

B   Ks

(1  Sr ) Sr pw 

B   Ks

(19)

(1  S r )2 pa

   kk       B (1  S r ) s  S r      ra  pa  a g     B (1  S r ) v  0 Ks    a 

(20)

where B1KT/Ks is the Biot’s coefficient, Ks and KT are the bulk moduli of soil particles and soil skeleton, respectively,  is the porosity, v is the volumetric strain, s is the suction, spapw, pw and pa are the pore water and pore air pressures (positive for compression), w and a are the densities of water and air, respectively, k is the intrinsic permeability tensor of soils, a and w are

10 / 46

ACCEPTED MANUSCRIPT

the dynamic viscosities of water and air, g is the gravitational acceleration vector, and Kw and Ka are the bulk moduli of water and air, respectively, with Kapa. As indicated in Section 2, since the water retention curve is deformation-dependent, the partial derivative of Sr with respect to time t in Eqs. (19) and (20) should be written as:

Sr 

S r S s r e s e

(21)

e  1  e0   v

CR IP T

The evolution of e is related to the volumetric strain by: (22)

Substituting Eqs. (21) and (22) into Eqs. (19) and (20) leads to:

AN US

 kk  cws v  cww pw  cwa pa      rw  pw   w g    0  w   kk  cas v  caw pw  caa pa      ra  pa  a g    0  a  with

(23)

(24)

cws   B S r 

M

deformation-dependent water retention

  B   S S r s  1  e0  r   Ks e  

ED

deformation-dependent water retention

PT

    S cas   B 1  S r    B 1  S r  s    1  e0  r e  Ks 

CE

cww 

AC

cwa  caw  caa 

B   Ks

B   Ks

B   Ks

B   Ks

S r2 

 Sr Kw



 S r  B   Sr s    s  Ks 

(25)

S r 1  S r  

 S r  B   Sr s    s  Ks 

S r 1  S r  

 S r  B   1  S r  s    s  Ks 

1  S r 

2



 1  S r  S r  Ka



  s 

B   Ks



1  S r  s  

It is worthy of noting that the coefficients cws and cas in Eqs. (23) and (24) are different from the previous studies because deformation-independent water retention curve models were always employed in these works [12, 51, 70]. However, as indicated in the coupled hydro-mechanical

11 / 46

ACCEPTED MANUSCRIPT

model, the deformation of soil skeleton affects the water retention behavior, which leads to an additional term (as labeled in Eq.(25)) in the two-phase flow governing equations. On the other hand, the governing equation for the solid deformation is given as:

     B  Sr pw  1  Sr  pa      g = 0

(26)

where  is the average skeleton stress tensor (positive for compression),  is the Kronecker delta

density of soil particles.

3.2 Finite element formation with modified coefficient matrix

CR IP T

tensor, and  is the total density of soils given by Srw+(1Sr)a+(1)s, with s being the

AN US

The displacement vector u(t)ux(t), uy(t), uz(t)}, pore water pressure pw(t) and pore air pressure pa(t) are selected as basic variables. By applying the finite element method into the governing equations (Eqs. (23), (24) and (26)), the discretize matrix equations for the coupled processes





BuT d  Csw pw  Csa pa  fs

(27)

M

can be obtained:

(28)

ED

Cws u  Cww pw  Cwa pa  K ww pw  fw Cas u  Caw pw  Caa pa  Kaa pa  fa

(29)

PT

where u, pw and pa are the vectors of nodal displacement, pore water pressure and pore air pressure, respectively, and Bu is the geometric matrix. The coefficient matrices are given in

CE

Appendix A.1.

In Eqs. (27)-(29), since the water retention curve model is deformation-dependent (as

AC

represented by cws and cas), the coefficient matrices of Cws and Cas are modified as:

C ws

deformation - dependent water retention       B   S r T  m T B d   N p  B Sr  S r s  1  e0    u    Ks e      

12 / 46

(30)

ACCEPTED MANUSCRIPT deformation - dependent water retention      B    S r  T T Cas   N p  B 1  S r    m Bu d 1  S r  s    1  e0    Ks e       

(31)

where m[1 1 1 0 0 0]T for three dimensional stress state and m[1 1 0]T for plain strain state, Np is the shape function.

CR IP T

The coupled matrix equations expressed by Eqs. (27)-(31) are highly nonlinear due to the strong nonlinearities of both constitutive relations and governing equations. This high nonlinearity would lead to complexity of the differentiation and more computational effort if a global Newton-Raphson procedure with consistent tangent moduli is employed. Here, following

AN US

the procedure reported by Schrefler and Scotta [51], a modified Newton-Raphson method using the elastic stiffness matrix rather than the tangent moduli matrix for predicting soil skeleton deformation is adopted to solve these coupled equations. By applying the generalized- integration scheme to discretize the coupled matrix equations in the time domain [71], the

M

incremental primary variables X = [u pw pa]T at the m-th iteration and n+1-th time step can be determined by solving the linearized coupled equations [51]:

ED

   m m m m m  t 1   f  Cn 1 + 1   f  K n 1  X n 1  rext,n 1  rint,n 1  

PT

with

AC

CE

C nm1

m rext, n 1

 0  C ws   Cas

0 C ww Caw

m

0   K ss  m C wa , K n 1   0    0 Caa  n 1

Csw K ww 0

(32)

m

Csa  0   Kaa  n 1

m  BuT nm1d   0   fs        m m   fw  , rint,  0  C nm1 X nm1  K nm1  pw,   n +1   n  1     m     0  fa  n 1  pa,n 1   

(33)

(34)

In Eqs. (32)-(34), m, f,  and  are the time integration parameters [72] and the matrix Kss is given by Kss   BuT D e Bu d , where De is the matrix form of the fourth-order elastic stiffness 

tensor Ce, known as the elastic moduli matrix. Unlike the the global Newton-Raphson procedure with consistent tangent moduli [49], the modified Newton-Rapson scheme does not possess a 13 / 46

ACCEPTED MANUSCRIPT

quadratic convergence rate. Nevertheless, a satisfactory convergence rate is observed in our numerical examples, as shown later in sections 5 and 6. The global residual force m m rext, n 1  rint,n 1

2

m and the incremental values of the primary variables X n 1

2

are used as the

convergence criteria to terminate the iterations: 2

and

X nm1 2   dof X nm11  X n

where T and dof are the error tolerances.

2

CR IP T

m m m rext, n 1  rint,n 1 2   T rext,n 1

(35)

Note that the coefficient matrices of Cws and Cas (Eqs. (32)-(34)) are different from the previous works. With only consideration of Sr = Sr(s), Hu et al. [29] derived these matrices as follows:

AN US

Cws   N pT B Sr mT Bu d 

Cas   N pT B 1  Sr  mT Bu d 

(36) (37)

The convergence property of the numerical formulation with the modified coefficient matrix

M

(Eqs. (30)-(31)) will be compared to the original coefficient matrices (Eqs. (36)-(37)) in Section

ED

5.3.

PT

4 Return mapping algorithm

The above-mentioned coupled hydro-mechanical model includes six stress variables and an

CE

additional scalar variable  for the mechanical component, as well as suction s and void ratio e for the hydraulic component. In the incremental finite element formulation, the pore water

AC

pressure pw and air pressure pa are primary variables and their increments (dpw, dpa) are known at the current m-th iteration and n-th time step, thus the current state of suction s and its increment (ds) are known. As expressed by Eqs. (4) and (5), the incremental Sr can only be determined when ds and de are both known. However, the incremental void ratio de is related to the elasto-plastic behavior, and this elasto-plastic response is not only related to the stress variable (), but also closely related to the water retention behavior in turn. This coupling of the mechanical and hydraulic behaviors challenges the constitutive integration procedure. To update 14 / 46

ACCEPTED MANUSCRIPT

the stress-like variables and the hydraulic variable (Sr), we develop a return mapping algorithm to integrate the hydraulic and mechanical relations with cross iteration. At the n-th time step and the i-th iteration, the stress tensor ni, the plastic strain tensor np, the yield stress pc(0)n, the void ratio en, and the degree of saturation Sr, n are known. The procedure of the integration is to update these variables from the n-th to n+1-th time step during

CR IP T

which the Kuhn-Tucker loading-unloading conditions must be satisfied for the mechanical model and the hydraulic path should be located in the bounding surfaces defined by Eq. (3). Therefore, we have

(38)

AN US

  n   f ( , , pc (0))0   n 1  d 0 p   p      n  d f ( , , pc (0)) 0   n 1        Se,w  Se  Se,d  pc (0)n   pc (0) n 1     Se,n    Se,n 1  

We develop a two-step procedure to determine the incremental values expressed by Eq.(38),

M

including dn, dnp, dpc(0)n, dSr, n.

4.1 Integration procedure for the hysteretic water retention behavior

ED

Suppose that the current hydraulic states are in the scanning zone, according to Eq.(5), the

PT

incremental saturation is given by

dSe   Se (1  Se1 mw )  ( kss

ds  ksede) s

(39)

CE

To obtain the analytical solution of Eq. (39), we arrange it with

AC

Se1 mw 1dSe  1 mw  kssd ln s  ksede Se 1

(40)

Eq. (40) can be further simplified into:

d ln  Se1 mw  1

mw

 kssd ln s  ksede

(41)

Integrating Eq. (41) from Se, n to Se, n+1, sn to sn+1, and en to en+1, yields

mw 

Se,n+1

Se,n

d ln  Se1 mw  1  kss  d ln s  kse  de sn1

en1

sn

en

(42)

Given that sn+1 = sn+dsn, en+1 = en+den and Se,n are known, the unknown variable Setrial, n+1

15 / 46

ACCEPTED MANUSCRIPT

can be calculated by:

   kss  kss mw  sn 1 exp  en 1   mw   1 mw  trial  Se,n 1  1  1  Se,n     k k m ss  snss w exp  en    mw   

 mw

(43)

CR IP T

The updated saturation determined by Eq. (43) is a trial saturation because it should be located in the bounding surfaces defined by Se,w(sn+1, en+1) and Se,d(sn+1, en+1) (Eq. (3)). The following procedure is proposed to correct the trial saturation:

trial trial trial If Se,n 1  Se,d ( sn 1 , en 1 ) and Se,n1  Se,w ( sn1 , en1 ) , then Se,n 1  Se,n 1 ;

else if Se,n 1  Se,d ( sn 1 , en 1 ) , then Se,n1  Se,d ( sn1 , en1 ) ;

AN US

trial

else Se,n1  Se,w ( sn1 , en1 ) , then Se,n1  Se,w ( sn1 , en1 ) . trial

4.2 Integration procedure for coupled hydro-mechanical constitutive model

M

As motivated by Eq. (38), the increments of the stress state variables can be determined by

ED

solving the following equations composed of flow rule, hardening rule and consistent condition

 np1   np   n 1

CE

PT

[e.g., 29, 49],

g 0   n 1

 1 e g  pc (0)n 1  pc (0)n exp   n 1 0 p n 1     f ( n 1 , n 1 , pc (0)n 1 )  0

(44)

AC

together with the elastic law to make the above equations closed:

 n 1  Ce :   n 1   np1 

(45)

By the general closest point projection [43, 73], the strain tensor at n+1-th time step, n+1,

remains constant, and from Eq. (45) we have:  np1  Ce1 :  n 1

(46)

The above-mentioned cross iteration is denoted by i, with i 0, 1, 2…. At the i-th iteration

16 / 46

ACCEPTED MANUSCRIPT

step, the solution of Eqs. (44)-(46) yields the increments of the stress state variables:

 g i      n 1   Ri    ni1  1  p,n +1  2 i 1        i     n 1   i i pc (0)n 1   Rh,n +1   g   p   n 1 

(47)

f (  2  in 1 

i n 1

,

i n 1

, pc (0)

 f    

i n 1

 f )    i

CR IP T

with f  pc (0) 

i

 Rp,i n +1      i   Rh,n +1    n 1

 g f  g   1    pc (0)  n 1 p    

1

T i

(48)

n 1

AN US

The increments of the stress state variables (Eq. (47)) are related to the current hydraulic variable of Sr, because Eq. (48) includes the bonding factor which is a function of Sr. In fact, before calculating the increments by Eq. (47), the elastic predictor must be performed to check if the current stress state is on the yield surface. Before conducting this elastic predictor checking,

M

we incorporate the integration procedure of the water retention model into the mechanical

ED

integration procedure. Therefore, the stress state is updated according to the water retention behavior, which represents the coupling of hydraulic behavior to the mechanical behavior.

PT

Finally, we present the integration procedure for the coupled hydro-mechanical constitutive

CE

model, as shown in Table 1.

AC

5 Validation

In this section, the developed numerical method is evaluated by two existing laboratory tests to demonstrate the efficiency in describing the coupled hydro-mechanical behaviors of unsaturated soils. The numerical method regarding the coefficient matrices (given in Section 3.2) is later evaluated by a computational case of rain infiltration into a soil column.

5.1 Comparison with laboratory tests: isotropic loading together with a drying-wetting cycle

17 / 46

ACCEPTED MANUSCRIPT

The laboratory test by Raveendiraraj [74] is used to validate the numerical method, with the stress paths shown in Fig. 3. Note that the data has been used to validate the coupled hydro-mechanical constitutive model in Hu et al. [22], and here we focus on the performance of the developed numerical methods with and without hydraulic hysteresis. Using a similar calibration procedure, the model parameters are presented in Table 2. In these two cases (with and without hysteresis), all

CR IP T

of the model parameters are the same except the parameters d, d, kss and kse. When hydraulic hysteresis is neglected, the hydraulic path is assumed to remain on the main wetting/drying surfaces, and hence the parameter d or d can be determined by the initial hydraulic state. From the experimental data, the initial hydraulic state is: Sr, i = 0.602, ei = 1.117 and si = 300 kPa, which leads to d (or w) = 1.3469106 kPa-1. In the numerical simulations, a number of 100

AN US

incremental loading steps are specified for each stage (AB, BC, CD, DE, EF), and a total of 500 time steps are calculated. The convergence criteria are taken as T = 1.0103 and dof = 1.0104.

Figs. 4 and 5 present comparisons between the experimental data and the simulated results

M

with and without hydraulic hysteresis. For the mechanical behavior, both the elastic responses

ED

(AA, BC, CD, DE and EE) and the elasto-plastic responses (AB and EF) can be simulated by the numerical methods with and without hydraulic hysteresis. However, for the

PT

hydraulic behavior, due to the nature of water retention curve hysteresis, the predicted hydraulic behavior on the wetting stage (CD) is irreversible, whereas it is reversible on the drying stage

CE

(DE). The predicted saturation at the point of D is much lower than the observed one (Fig. 5d). This is because at D the soil sample underwent an equalization stage during which the suction

AC

and the net stress remained constant while the saturation changed [74], leading to a remarkable increase in saturation Sr and a slight increase in void ratio e. The hydraulic behavior in the equalization stage cannot be predicted by our constitutive hydromechanical model [22, 59] because the relationship for the elastic behavior (Eq. (13)) is not adequate for the expansive soil [59]. Nevertheless, as shown in Fig. 5d, the numerical method with hydraulic hysteresis can reasonably capture this hysteretic behavior, and the simulated results agree better with the measurements (Figs. 5c, d).

18 / 46

ACCEPTED MANUSCRIPT

Figs. 5a, b also show that hydraulic hysteresis occurring on the wetting-drying cycle (CDE) does not significantly affect the stress-strain behavior. The reason is that during this cycle the mechanical behavior is elastic, and the soil skeleton deformation is only governed by the elastic law (Eq. (12)). The change in degree of saturation during this cycle only slightly modifies the average skeleton stress and thus the variations of void ratio predicted by the two

CR IP T

numerical models are nearly the same. Nevertheless, hydraulic hysteresis significantly impacts the stress path because the bonding factor is sensitive to the degree of saturation, as shown in Fig. 6. For the numerical method with hydraulic hysteresis, a higher degree of saturation at point E is predicted, which induces a lower bonding factor. Thus, after point E, a remarkable elasto-plastic deformation is produced even though the mean net stress is smaller than its maximum value

AN US

previously applied (point E in Fig. 4b), which is an important behavior for unsaturated soils [17, 59]. Fig. 7 presents the convergence profiles of the proposed numerical method with hydraulic hysteresis at points B, E’ and F. The global residual force (T = 1.0103) converged within 25 global iterations, while the tolerances for the consistent condition (f1.0104) and the

M

flow/hardening rules (ph1.0106) could be achieved within 4 local iterations, indicating an

ED

overall good convergence property of the proposed numerical method.

PT

5.2 Comparison with laboratory tests: isotropic loading together with shearing The laboratory test performed by Raveendiraraj [74] including isotropic loading and shearing is

CE

employed to further validate the proposed methods. In the experiment, the soil sample was subjected to isotropic loading under a constant suction of 300 kPa as the mean net stress pnet

AC

increased from 10 to 75 kPa (AB), followed by subsequent shearing (BH) under fully drained conditions (q/pnet = 3) until a critical state was approached. All of the model parameters are the same as those given in Section 5.1, as presented in Table 2, except the model parameter d (or w). For the case of numerical simulation without hydraulic hysteresis, the model parameter d (or d) is only related to the initial state. From the experiment, the initial degree of saturation and the initial void ratio are respectively 0.606 and 1.15, which provides d (or w) = 6.0875  107 kPa-1. In the numerical simulations, a number of 100 incremental loading

19 / 46

ACCEPTED MANUSCRIPT

steps are specified for each stage and a total of 200 time steps are calculated. Fig. 8 compares the experimental data and the simulated results obtained by the numerical methods with and without hydraulic hysteresis, under different convergence criteria, with the global residual force error εT being specified in 10-2, 10-3 and 10-4, respectively. It illustrates that the predicted results obtained by the proposed method with hydraulic hysteresis agree better with

CR IP T

the measurements. Before A’, the mechanical behavior is elastic, and the predicted mechanical (Figs. 8, f) and hydraulic responses (Figs. 8c, d) for the numerical methods with and without hydraulic hysteresis are nearly the same. However, after A’, as the elasto-plastic deformation continues to occur, the degree of saturation is overestimated (Fig. 8c) by the numerical method without hydraulic hysteresis whereas a better match can be obtained when considering hydraulic

AN US

hysteresis (Fig. 8d). As the degree of saturation continues to be overestimated (Fig. 8c) in the case of neglecting hydraulic hysteresis, the computed stress path (Fig. 9) evolves much more sharply, and thus the expanded yield surface (C1) is larger. The convergence profile for the proposed numerical method with hydraulic hysteresis at the time steps (#75, #95, and #125) are

M

given in Fig. 10. The curves show that the convergence for the global residual force (εT) and for

ED

the local residuals could be achieved within 55 iterations and 4 iterations, respectively,

PT

indicating an overall good convergence property of the proposed numerical method.

5.3 Rain infiltration into a soil column

CE

In this section, the coupled flow-deformation process induced by rain infiltration into a soil column is simulated to investigate the performance of the numerical formation with modified

AC

coefficient matrices shown in Section 3.2. The soil column is 1 m in height and 0.1 m in diameter (Fig. 11). The initial effective degree of saturation is assumed to decrease linearly from 1.0 at z 0 to 0.6 at the top of the soil column (z  1). The initial pore air pressure is prescribed as the atmosphere pressure. The lateral surface is impermeable for water/air flow and fixed with zero normal displacement. The initial stress field is obtained by solving the equilibrium equation under gravity and the above initial/boundary conditions. A pre-load of 15 kPa is applied on the top surface to make the soil a little over-consolidated and to obtain the initial yield stress pc(0) at

20 / 46

ACCEPTED MANUSCRIPT

each Gauss point. Afterwards (t > 0), the pre-load is removed and the nodal displacements are set back to zero. A rainfall with intensity of 14.4 mm/h and a surface load of 10 kPa are then applied on the top surface (t > 0). The initial and maximum time step lengths are set to 0.001 h and 0.1 h, respectively, with an incremental factor of 1.05 under the convergence conditions of

T1.0E-3, dof1.0E-3,f1.0E-3 and ph1.0E-5. The material properties are listed in Table 3, and the constitutive model parameters are given in Table 4.

CR IP T

Figs. 12a, b and c show the simulated results at z  0.8 m predicted by the numerical methods with the original and the modified coefficient matrices. Before point A, the mechanical response is elastic and the numerical results with the original and modified coefficient matrices are nearly the same. However, after point A, a fluctuation of air pressure is observed (Fig. 12b)

AN US

when using the original Cws and Cas, and a remarkable decline of the vertical displacement was predicted in this case (Fig. 12c). As for the simulated results at z  0.2 m (Figs. 12d, e and f), the evolutions of water pressure for the two cases are nearly the same, due to a longer distance to the water infiltration boundary (top surface). Again, a fluctuation of air pressure is found at t = 80

M

min.

ED

Figs 12a~f also present the simulated results with the three different meshes (Figure 11). As the modified coefficient matrices were employed, the evolutions of the primary variables, i.e. pw,

PT

pa and uz, are nearly the same for the coarse, medium and fine meshes. Table 5 summarizes the simulated results of pw, pa and uz at the end of simulations at z = 0.2 m and 0.8 m, which shows

CE

that the relative differences of pw, pa and uz to the fine mesh are no more than 3.02%. Thus, the mesh refinement shown in Figure 11 has no significant influence on the results, which indicates

AC

that the mesh-insensitive results can be obtained by the numerical method with the modified Cws and Cas using the coarse mesh. Figs. 12g and h show the convergence profile for the global residual force (T) at the time

steps of #75and #76, corresponding to t = 69.3 min and t = 74.1 min, respectively. It indicates that for the two cases the global convergence with εT = 1.0×10-3 can both be achieved within 20 global iterations at time steps of #75. But at the time step of #76, Fig. 12h shows a better convergence rate of the numerical method with the modified Cws and Cas. Only 18 global iterations are needed to

21 / 46

ACCEPTED MANUSCRIPT

obtain the convergent results, instead of 80 global iterations in the case using the original Cws and Cas.

6 Application

CR IP T

In this section, the proposed numerical method is applied for analysis of the coupled water-air two-phase flow and elasto-plastic deformation processes in a homogenous soil slope under rain infiltration, aiming to gain a better understanding of the role of hydraulic hysteresis in the coupled water-air flow and deformation processes.

AN US

6.1 Computational model

Fig. 13 shows a soil slope model of 10 m high, with a gradient of 2:3 and an initial groundwater level of 0 m. The finite element meshes of the computational domain include 1434 four-node quadrilateral elements and 1514 nodes. All of the parameters included in the governing equations

M

are listed in Table 3, and the parameters for the coupled hydro-mechanical constitutive model are

ED

presented in Table 6. As mentioned previously, the main drying curve for water retention behavior is always used for the coupled hydro-mechanical analysis of the rainfall-induced slope deformation. Thus for comparison, the model parameter d = w is specified as 0.02857 kPa-1 by

PT

neglecting the nature of hydraulic hysteresis. In the case of considering hydraulic hysteresis, d =

CE

0.02857 kPa1, w = 0.07931 kPa1, kss = 0.38 and kse = 0.58. The initial effective degree of saturation is assumed to decrease linearly from 1.0 at the initial groundwater level (z = 0) to 0.6 at

AC

the crest of the slope (z = 0.15). The initial pore air pressure (pa) is prescribed as the atmospheric pressure. The boundary BC is specified as the Dirichlet boundary with pw=0 and pa=patm, and fixed with zero horizontal and vertical displacement; the boundaries of AB and CD are assumed to be impervious for water/air flow, and fixed with zero horizontal displacement; the boundaries of AF, FE and ED are subjected to rain infiltration for water flow and specified with pa=patm. The initial stress field is obtained by the solution of the equilibrium equation under gravity and the above initial/boundary conditions. A preload of 10 kPa is then applied on the slope surface (AF,

22 / 46

ACCEPTED MANUSCRIPT

FE and ED) to make the soil slightly over-consolidated and consequently determine the initial yield stress pc(0) at each Gauss point. Afterward, the preload and the nodal displacements are removed. The slope surface is subjected to a rain infiltration of 16 mm/h. The initial and maximum time step lengths are set to 0.001 and 0.05 h, respectively, with the following

CR IP T

convergence conditions: T1.0E-3, dof1.0E-3,f1.0E-3 and ph1.0E-5.

6.2 Computational results

Fig. 14 shows the distributions of the plastic deviatoric strain (qp) calculated by the water retention curve model with (Figs. 14a, c and e) and without hydraulic hysteresis (Figs. 14 b, d and f) at t = 6, 7 and 7.3 h. Figs. 14a-f indicate that a shear band where the shear strain is

AN US

localized over a narrow zone forms at the toe of the slope and then propagates toward the slope crest, which agrees with the experimental observations [32]. Comparisons of the calculated shear bands presented in Figs. 14a, c and e and in Figs. 14b, d and f demonstrate that hydraulic hysteresis of water retention curve has a significant effect on the evolution of the shear band

M

during rain infiltration: the shear band calculated by the non-hysteretic water retention curve

ED

develops “faster” than that by the hysteretic water retention curve. The simulated results of qp at the failure state are shown in Figs. 14g, h, in which the slope failure is regarded to occur as the

PT

iteration does not converge (see Figure 17). The slope failure occurs at 7.5 h when hydraulic hysteresis is neglected. But the time of failure would delay by 0.9 h if hydraulic hysteresis is

CE

taken into account.

Figs. 15 and 16 further indicate the considerable impact of hydraulic hysteresis on the

AC

coupled flow-deformation processes. At the failure state of the slope, the plastic deviatoric strain is accumulated within a narrow zone (Figs. 14g, h), which affects the distributions of the displacement vector and the degree of saturation. For the simulated results without consideration of hydraulic hysteresis, the maximum displacement occurs near the center of the slope surface (EF), with a magnitude of 0.63 m (Fig. 15a), whereas in the second case (with hydraulic hysteresis) the maximum displacement is 0.27 m and occurs at the toe of slope (Fig. 15b). On the other hand, given the deformation-dependent water retention curve considered in the simulation, the

23 / 46

ACCEPTED MANUSCRIPT

distribution of the degree of saturation also exhibits a zoning feature. Figs. 15c, d show a higher degree of saturation in the zone, induced by the decrease of void ratio. Because larger deformation is predicted by the non-hysteretic water retention curve, a larger zone of the degree of saturation is produced. The impact of hydraulic hysteresis on the coupled flow-deformation can also be illustrated

CR IP T

by the evolutions of the horizontal displacement, vertical displacement, degree of saturation and the plastic deviatoric strain at a typical point P (see its location in Fig. 13), as shown in Fig. 16. It can be observed that the predicted evolutions of these variables are all delayed with lower values by the hysteretic water retention curve. For the case of non-hysteretic water retention curve, the mechanical response is elastic before t = 6.88 h (A1) and after A1 plastic deformation

AN US

is produced, whereas the plastic deformation would occur at t = 7.83 (A2) in the case of hysteretic water retention curve. Figure 17 shows the convergence profiles of the proposed numerical method with hydraulic hysteresis at t = 3 h, 6 h and 8.4 h, respectively. At t = 3 h and 6 h, the specified error tolerances (T1.0E-3, dof1.0E-3) could be achieved within 80 global

M

iterations, but at the time of 8.4 h, the global residual force and the incremental value of the

ED

primary variables cannot converge, indicating the occurrence of the slope failure. Fig. 18 shows the predicted stress path and hydraulic path during the rain infiltration. It is

PT

demonstrated that the soil at point A experiences a softening process because the stress path during the whole rain infiltration process is located on the left side of the point at which the

CE

critical state line intersects the yield surface (known as dry side). Before the yielding point (A1 and A2), relatively small deformation and saturation are produced. After these points, substantial

AC

deformation is produced and the degree of saturation increases considerably. However, in the case of non-hydraulic hysteresis, the compressive deformation and the degree of saturation would be overestimated. Consequently, the yield stress at the final state (failure) in the case of non-hysteretic water retention curve is 51.2 kPa whereas this value is 66.1 kPa in the case of consideration of hydraulic hysteresis. The final yield surfaces for the two cases are also shown in Fig. 18. Note that even previous investigations showed that the occurrence of shallow landslides (represented by the factor of safety) would be underestimated using the main drying curve, our

24 / 46

ACCEPTED MANUSCRIPT

numerical results imply that the predicted shear band and the slope deformation would be overestimated by the non-hysteretic water retention curve. Given that the factor of safety is a simplified indicator without consideration of the solid-liquid coupling, it is preferable to evaluate the slope stability in the framework of coupled two-phase fluid flow and elasto-plastic

CR IP T

deformation.

7 Conclusions

We proposed a numerical method to examine the impact of hydraulic hysteresis of water retention curve on the coupled water-air two-phase flow and elasto-plastic deformation processes.

AN US

An implicit scheme was developed for integrating the stress-strain model and the hysteretic water retention curve model. We modified the coupled coefficient matrices to account for the deformation-dependent water retention curve. The proposed numerical formulation was first validated by two existing laboratory tests, indicating that the proposed method with

M

consideration of hydraulic hysteresis is able to better capture the coupled hydro-mechanical

ED

behavior under the various loading conditions, including isotropic loading, drying-wetting cycle and shearing. The computational example of rain infiltration into a soil column shows that the

PT

proposed numerical method with modified coefficient matrices has a better convergence property than that with the original coefficient matrices. Finally, we employed the proposed

CE

numerical method to examine the impacts of hydraulic hysteresis on the coupled water-air flow and deformation processes in a soil slope under rain infiltration. The numerical results show that

AC

the proposed procedure not only captures the development and evolution of shear band, but also reveals the significant effects of hydraulic hysteresis on the flow and deformation processes. The neglect of hydraulic hysteresis of the water retention curve would significantly overestimate the slope deformation, with an early prediction of the occurrence of failure.

Acknowledgments The financial supports from the National Natural Science Foundation of China (Nos. 51409198, 51579188) and China Postdoctoral Science Special Foundation (No. 2015T80833) are gratefully 25 / 46

ACCEPTED MANUSCRIPT acknowledged.

APPENDIX A A.1 The coefficient matrix involved in the governing equation (Eqs. (27)-(29)).

   T Csa   Bu  1  S r  mN pd    fs   N uT  gd   N uT t d   

(A1)

  T  C wa   cwa N p N pd    k k  K ww   (N p )T rw (N p )d   w  T kk rw T fw   (N p ) gd   q N p qw d   w  w

(A2)

  T  Caa   caa N p N pd    T kk ra K aa   (N p ) (N p )d    a  T kk ra T fa   (N p ) gd   q N p qa d   a  a

(A3)

Csw   BuT S r mN pd

C ww   cww N pT N pd

AN US



CR IP T



Caw   caw N pT N pd

R.I. Borja, J.A. White, X. Liu, W. Wu, Factor of safety in a partially saturated slope inferred from hydro-mechanical continuum modeling, Int. J. Numer. Anal. Methods Geomech., 36 (2012) 236-248. S.E. Cho, S.R. Lee, Instability of unsaturated soil slopes due to infiltration, Comput. Geotech., 28 (2001) 185-208. L. Laloui, A. Ferrari, C. Li, J. Eichenberger, Hydro-mechanical analysis of volcanic ash slopes during rainfall, Geotechnique, 66 (2016) 220-231. R.I. Borja, X. Liu, J.A. White, Multiphysics hillslope processes triggering landslides, Acta Geotech., 7 (2012) 261-269. E. Liu, H.-S. Yu, G. Deng, J. Zhang, S. He, Numerical analysis of seepage–deformation in unsaturated soils, Acta Geotech., 9 (2014) 1045-1058. M. Raj, A. Sengupta, Rain-triggered slope failure of the railway embankment at Malda, India,

AC

1.

CE

References

PT

ED

M



2. 3. 4. 5. 6.

26 / 46

ACCEPTED MANUSCRIPT

13. 14. 15.

16.

17. 18.

AC

19.

CR IP T

12.

AN US

11.

M

10.

ED

9.

PT

8.

CE

7.

Acta Geotech., 9 (2014) 789-798. A. Lloret, M. Villar, Advances on the knowledge of the thermo-hydro-mechanical behaviour of heavily compacted “FEBEX” bentonite, Physics and Chemistry of the Earth, Parts A/B/C, 32 (2007) 701-715. A. Gens, M. Sánchez, L.D.N. Guimaraes, E. Alonso, A. Lloret, S. Olivella, M. Villar, F. Huertas, A full-scale in situ heating test for high-level nuclear waste disposal: observations, analysis and interpretation, Geotechnique, 59 (2009) 377-399. J. Rutqvist, Y. Ijiri, H. Yamamoto, Implementation of the Barcelona Basic Model into TOUGH–FLAC for simulations of the geomechanical behavior of unsaturated soils, Comput. Geosci., 37 (2011) 751-762. M. Sánchez, A. Gens, L. Guimarães, Thermal–hydraulic–mechanical (THM) behaviour of a large-scale in situ heating experiment during cooling and dismantling, Can. Geotech. J., 49 (2012) 1169-1195. M.M.d. Farias, M.P. Cordão Neto, Advanced numerical simulation of collapsible earth dams, Can. Geotech. J., 47 (2010) 1351-1364. A.R. Khoei, T. Mohammadnejad, Numerical modeling of multiphase fluid flow in deforming porous media: A comparison between two- and three-phase models for seismic analysis of earth and rockfill dams, Comput. Geotech., 38 (2011) 142-166. A. Khoei, E. Haghighat, Extended finite element modeling of deformable porous media with arbitrary interfaces, Appl. Math. Model., 35 (2011) 5426-5441. E.E. Alonso, A. Gens, A. Josa, A constitutive model for partially saturated soils, Geotechnique, 40 (1990) 405-430. J. Vaunat, E. Romero, C. Jommi, An elastoplastic hydromechanical model for unsaturated soils, in: Trento (Ed.) Experimental evidence and theoretical approaches in unsaturated soils, 2000, pp. 121–138. D. Gallipoli, A. Gens, R. Sharma, J. Vaunat, An elasto-plastic model for unsaturated soil incorporating the effects of suction and degree of saturation on mechanical behaviour, Geotechnique, 53 (2003) 123-136. S. Wheeler, R. Sharma, M. Buisson, Coupling of hydraulic hysteresis and stress-strain behaviour in unsaturated soils, Geotechnique, 53 (2003) 41-54. R. Tamagnini, An extended cam-clay model for unsaturated soils with hydraulic hysteresis, Geotechnique, 54 (2004) 223-228. R. Santagiuliana, B. Schrefler, Enhancing the Bolzon–Schrefler–Zienkiewicz constitutive model for partially saturated soil, Transport Porous. Med., 65 (2006) 1-30. D.a. Sun, D. Sheng, S.W. Sloan, Elastoplastic modelling of hydraulic and stress–strain behaviour of unsaturated soils, Mech. Mater., 39 (2007) 212-221. D. Sheng, D.G. Fredlund, A. Gens, A new modelling approach for unsaturated soils using independent stress variables, Can. Geotech. J., 45 (2008) 511-534. R. Hu, Y.-F. Chen, H.-H. Liu, C.-B. Zhou, A coupled stress–strain and hydraulic hysteresis model for unsaturated soils: Thermodynamic analysis and model evaluation, Comput. Geotech., 63 (2015) 159-170. A. Zhou, D. Sheng, An advanced hydro-mechanical constitutive model for unsaturated soils with different initial densities, Comput. Geotech., 63 (2015) 46-66.

20. 21. 22.

23.

27 / 46

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

24. M. Nuth, L. Laloui, Effective stress concept in unsaturated soils: clarification and validation of a unified framework, Int. J. Numer. Anal. Methods Geomech., 32 (2008) 771-801. 25. K.K. Muraleetharan, C. Liu, C. Wei, T.C.G. Kibbey, L. Chen, An elastoplatic framework for coupling hydraulic and mechanical behavior of unsaturated soils, Int. J. Plasticity., 25 (2009) 473-490. 26. A.-N. Zhou, D. Sheng, S.W. Sloan, A. Gens, Interpretation of unsaturated soil behaviour in the stress–Saturation space, I: Volume change and water retention behaviour, Comput. Geotech., 43 (2012) 178-187. 27. M. Nuth, L. Laloui, Advances in modelling hysteretic water retention curve in deformable soils, Comput. Geotech., 35 (2008) 835-844. 28. B. Loret, N. Khalili, A three-phase model for unsaturated soils, Int. J. Numer. Anal. Methods Geomech., 24 (2000) 893-927. 29. R. Hu, Y.-F. Chen, H.-H. Liu, C.-B. Zhou, A coupled two-phase fluid flow and elastoplastic deformation model for unsaturated soils: theory, implementation, and application, Int. J. Numer. Anal. Methods Geomech., 40 (2016) 1023-1058. 30. D. Sheng, S.W. Sloan, A. Gens, D.W. Smith, Finite element formulation and algorithms for unsaturated soils. Part I: Theory, Int. J. Numer. Anal. Methods Geomech., 27 (2003) 745-765. 31. F. Oka, S. Kimoto, N. Takada, H. Gotoh, Y. Higo, A seepage-deformation coupled analysis of an unsaturated river embankment using a multiphase elasto-viscoplastic theory, Soils Found, 50 (2010) 483-494. 32. Y. Xiong, X. Bao, B. Ye, F. Zhang, Soil–water–air fully coupling finite element analysis of slope failure in unsaturated ground, Soils Found, 54 (2014) 377-395. 33. A. Gens, Soil-environment interactions in geotechnical engineering, Geotechnique, 60 (2010) 3-74. 34. D. Sheng, A. Gens, D.G. Fredlund, S.W. Sloan, Unsaturated soils: from constitutive modelling to numerical algorithms, Comput. Geotech., 35 (2008) 810-824. 35. D. Sheng, S. Sloan, A. Gens, A constitutive model for unsaturated soils: thermomechanical and computational aspects, Comput. Mech., 33 (2004) 453-465. 36. F. Nagel, G. Meschke, An elasto-plastic three phase model for partially saturated soil for the finite element simulation of compressed air support in tunnelling, Int. J. Numer. Anal. Methods Geomech., 34 (2010) 605-625. 37. N. Khalili, M.A. Habte, S. Zargarbashi, A fully coupled flow deformation model for cyclic analysis of unsaturated soils including hydraulic and mechanical hystereses, Comput. Geotech., 35 (2008) 872-889. 38. J. Vaunat, J. Cante, A. Ledesma, A. Gens, A stress point algorithm for an elastoplastic model in unsaturated soils, Int. J. Plasticity., 16 (2000) 121-141. 39. W. Ehlers, T. Graf, M. Ammann, Deformation and localization analysis of partially saturated soil, Comput. Method. Appl. M, 193 (2004) 2885-2910. 40. A. Khoshghalb, N. Khalili, A meshfree method for fully coupled analysis of flow and deformation in unsaturated porous media, Int. J. Numer. Anal. Methods Geomech., 37 (2013) 716-743. 41. X. Song, R.I. Borja, Mathematical framework for unsaturated flow in the finite deformation range, Int. J. Numer. Meth. Eng., 97 (2014) 658-682. 28 / 46

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

42. B. Shahbodagh-Khan, N. Khalili, G.A. Esgandani, A numerical model for nonlinear large deformation dynamic analysis of unsaturated porous media including hydraulic hysteresis, Comput. Geotech., 69 (2015) 411-423. 43. J.C. Simo, R.L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Comput. Method. Appl. M, 48 (1985) 101-118. 44. S. Sloan, Substepping schemes for the numerical integration of elastoplastic stress–strain relations, Int. J. Numer. Meth. Eng., 24 (1987) 893-911. 45. J. Li, S. Chen, L. Jiang, On implicit integration of the bounding surface model based on swell–shrink rules, Appl. Math. Model., 40 (2016) 8671-8684. 46. T. Ma, C. Wei, H. Wei, W. Li, Hydraulic and Mechanical Behavior of Unsaturated Silt: Experimental and Theoretical Characterization, Int. J. Geomech., 16 (2016) D4015007. 47. T.-T. Ma, C.-F. Wei, P. Chen, H.-Z. Wei, Implicit scheme for integrating constitutive model of unsaturated soils with coupling hydraulic and mechanical behavior, Appl. Math. Mech., 35 (2014) 1129-1154. 48. Y. Zhang, A. Zhou, Explicit integration of a porosity-dependent hydro-mechanical model for unsaturated soils, Int. J. Numer. Anal. Methods Geomech., 40 (2016) 2353-2382. 49. H.W. Zhang, L. Zhou, Implicit integration of a chemo-plastic constitutive model for partially saturated soils, Int. J. Numer. Anal. Methods Geomech., 32 (2008) 1715-1735. 50. R.I. Borja, J.A. White, Continuum deformation and stability analyses of a steep hillside slope under rainfall infiltration, Acta Geotech., 5 (2010) 1-14. 51. B.A. Schrefler, R. Scotta, A fully coupled dynamic model for two-phase fluid flow in deformable porous media, Comput. Method. Appl. M, 190 (2001) 3223-3246. 52. T.-L. Tsai, Influences of soil water characteristic curve on rainfall-induced shallow landslides, Environ. Earth Sci., 64 (2010) 449-459. 53. K.-C. Ma, Y.-C. Tan, C.-H. Chen, The influence of water retention curve hysteresis on the stability of unsaturated soil slopes, Hydrol. Process., 25 (2011) 3563-3574. 54. S.E. Cho, Stability analysis of unsaturated soil slopes considering water-air flow caused by rainfall infiltration, Eng. Geol., 211 (2016) 184-197. 55. K.-H. Yang, R. Uzuoka, J.N.a.a. Thuo, G.-L. Lin, Y. Nakai, Coupled hydro-mechanical analysis of two unstable unsaturated slopes subject to rainfall infiltration, Eng. Geol., (2016). 56. R. Hu, Y. Chen, C. Zhou, Modeling of coupled deformation, water flow and gas transport in soil slopes subjected to rain infiltration, Sci. China Tech. Sci., 54 (2011) 2561-2575. 57. B.A. Ebel, K. Loague, R.I. Borja, The impacts of hysteresis on variably saturated hydrologic response and slope failure, Environ. Earth Sci., 61 (2010) 1215-1225. 58. R. Hu, Y.F. Chen, H.H. Liu, C.B. Zhou, A water retention curve and unsaturated hydraulic conductivity model for deformable soils: consideration of the change in pore size distribution., Geotechnique, 63 (2013) 1389-1405. 59. R. Hu, H.-H. Liu, Y. Chen, C. Zhou, D. Gallipoli, A constitutive model for unsaturated soils with consideration of inter-particle bonding, Comput. Geotech., 59 (2014) 127-144. 60. A. Tarantino, S. Tombolato, Coupling of hydraulic and mechanical behaviour in unsaturated compacted clay, Geotechnique, 55 (2005) 307-317. 61. A. Tarantino, A water retention model for deformable soils, Geotechnique, 59 (2009) 751-762. 29 / 46

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

62. R. Lewis, B. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. Wiley,Chichester, 1998. 63. A. Bishop, The principle of effective stress, Teknisk Ukeblad, 106 (1959) 859-863. 64. G. Houlsby, The work input to an unsaturated granular material, Geotechnique, 47 (1997) 193-196. 65. S.M. Hassanizadeh, W.G. Gray, Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries, Adv. Water Resour., 13 (1990) 169-186. 66. N. Khalili, M. Khabbaz, A unique relationship of chi for the determination of the shear strength of unsaturated soils, Geotechnique, 48 (1998) 681-687. 67. N. Khalili, F. Geiser, G. Blight, Effective stress in unsaturated soils: review with new evidence, Int. J. Geomech., 4 (2004) 115-126. 68. K.H. Roscoe, J. Burland, On the generalized stress-strain behaviour of wet clay, (1968). 69. J. Jaky, Pressure in silos, 1948, 103-107. 70. R. Hu, Y. Chen, H. Liu, C. Zhou, A relative permeability model for deformable soils and its impact on coupled unsaturated flow and elasto-plastic deformation processes, Sci. China Tech. Sci., 58 (2015) 1971-1982. 71. J. Chung, G. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method, J. Appl. Mech., 60 (1993) 371-375. 72. D. Kuhl, M. Crisfield, Energy-conserving and decaying algorithms in non-linear structural dynamics, Int. J. Numer. Meth. Eng., 45 (1999) 569-599. 73. J. Simo, T. Hughes, Computational inelasticity, Springer: New York1998. 74. A. Raveendiraraj, Coupling of mechanical behaviour and water retention behaviour in unsaturated soils, University of Glasgow, Ph.D. thesis, 2009.

30 / 46

ACCEPTED MANUSCRIPT

Figures (18 in total)

c

AN US

b

lnSe

DS

kss 1 kss 1

Constant e

1

DS

kse 1 kse 1 mnkp

mnkp

1

1

1

e

ED

lns

DM

Se =1

Constant s

mn

mn

M

lnSe

Se =1

DM

CR IP T

a

Fig. 1. The hysteretic water retention curve for deformable soils in the space of s-e-Se. (a) Main

AC

CE

PT

drying and wetting surfaces; scanning curve in Selns plane (b) and in See plane (c). q Initial yield surface M

pc'(0) 

p'

Subsequent yield surface

Fig. 2. The initial and subsequent yield surfaces of the elastoplastic model in the p--q space.

31 / 46

ACCEPTED MANUSCRIPT

A

F

B

C E

Stress path: ABCDEF

D 10 O 10

200 pnet (kPa)

CR IP T

s (kPa)

300

375

Fig. 3. Experimental stress path involving isotropic loading and a drying–wetting cycle.

a 1.25

D

D

1.15 A

1.00 C,E

D

0.70 C,E

0.65

10

B

E'

Measured Predicted elastic response Predicted elastoplastic response

d D F

Measured Predicted reversible response Predicted irreversible response

F

F

E B C

B A''

A

A

AC

0.60

F

Measured Predicted reversible response

CE

Sr (-)

0.75

ED

0.85 0.85

A'

C,E

B

Measured Predicted elastic response Predicted elastoplastic response

0.90

0.80

E'

M

0.95

c

A

A'

1.05

PT

e (-)

1.10

with hydraulic hysteresis

AN US

1.20

b

without hydraulic hysteresis

100

10

100

pnet (kPa)

pnet (kPa)

Fig. 4. Comparisons between experimental data (dot) and the predictions (line) by the coupled model without hydraulic hysteresis (a, c) and with hydraulic hysteresis (b, d) for the isotripic loading stage (ABC, EF) shown in Figure 3. (a, b) the change in the void ratio with the mean net stress; (c, d) the change in the degree of saturation with the mean net stress.

32 / 46

ACCEPTED MANUSCRIPT

a

b

without hydraulic hysteresis 1.25 1.20 D

with hydraulic hysteresis

D

1.15

e (-)

1.10 1.05 Measured Predicted elastic response

0.95

Measured Predicted elastic response

CR IP T

1.00

E C

E C

0.90

c 0.80

d

D D

0.70

E

E C

AN US

Sr (-)

0.75

0.65

Measured Predicted reversible response Predicted irreversible response

Measured Predicted reversible response

0.60 10

100

10

s (kPa)

C

100

s (kPa)

M

Fig. 5. Comparisons between experimental data (dot) and the predictions (line) by the coupled

ED

model without hydraulic hysteresis (a, c) and with hydraulic hysteresis (b, d) for the drying-wetting cycle (CD, DE) shown in Figure 3. (a, b) the change in the void ratio with the

0.024

AC

Predicted elastic response 0.022 Predicted 0.020 elastoplastic response

 (-)

b

without hydraulic hysteresis

CE

a

PT

suction; (c, d) the change in the degree of saturation with the suction.

A

E'

0.018

Predicted elastic response Predicted elastoplastic response

A'

C,E

0.016

with hydraulic hysteresis

B

A

A'

C

E F

B E' F

0.014 0.012

0.010 10

D D

C0 100

C1

C2 1000 10

p' (kPa)

C0 100

C1

C2 1000

p' (kPa)

Fig. 6. The computed stress paths involving isotropic loading and a drying–wetting cycle: (a)

33 / 46

ACCEPTED MANUSCRIPT

without hydraulic hysteresis, (b) with hydraulic hysteresis

101 B E' F

100 10-1 10-2 10-3 10-4 10-5

0

5

10 15 20 25 30 35 40 45 50

m, number of global iteraions

c

102

B E' F

100

10-4 10-6 10-8 10-10

B E' F

10-4 10-5 10-6 10-7 10-8 10-9

10-10

-12

1

2

3

M

10

10-3 i i Rp，  Rh，  n+1 1 n+1

i fn+1 (kPa2)

10-2

10-2

AN US

b

CR IP T

rm rm  /rm  ext， n+1 int， n+1 2 ext， n+1 2

a

4

10-11

2

3

i, number of local iteraions

ED

i, number of local iteraions

1

Fig. 7. Convergence profiles of the proposed numerical method with hydraulic hysteresis for the

PT

isotropic loading and shearing test at points B, E’ and F which are shown in Figures 3 and 4. (a) the global residual force, (b) the residual of consistent condition, and (c) the residual for the

AC

CE

hardening rule and flow rule.

34 / 46

ACCEPTED MANUSCRIPT

a

b

without hydraulic hysteresis

300

with hydraulic hysteresis H

H

250

150 Measured Predicted (=1.0E4)

50

Predicted (=1.0E3)

0.05

0.10

0.85

0.20

0.25

0.30 0.00

A'

A'

1.00

f

A'

H

B

Predicted (=1.0E3)

CE 10

B

Measured Predicted (=1.0E4) Predicted (=1.0E3)

H

Predicted (=1.0E2)

100

pnet (kPa)

A'

B

PT

0.90

0.30

A

Measured Predicted (=1.0E4)

0.95

100

10

pnet (kPa)

H

Predicted (=1.0E2)

100 pnet (kPa)

AC

0.80

10

ED

A

1.05

0.85

100

pnet (kPa)

e 1.20 1.10

0.25

A

A 10

1.15

AN US

0.60

0.20

Predicted (=1.0E2)

B

0.65

0.15

Predicted (=1.0E3)

Predicted (=1.0E2)

0.70

0.10

Measured Predicted (=1.0E4)

Predicted (=1.0E3)

0.75

0.05

q (-)

d

H

Measured Predicted (=1.0E4)

0.80 Sr (-)

0.15

Predicted (=1.0E2)

B

q (-)

c

e (-)

Predicted (=1.0E3)

Predicted (=1.0E2)

B

0 0.00

Measured Predicted (=1.0E4)

CR IP T

100

M

q (kPa)

200

HFig. 8. Comparisons between experimental data (dot) and predictions (line) by the coupled model

without hydraulic hysteresis (a, c, e) and with hydraulic hysteresis (b, d, f) for the stress path involving isotropic loading and shearing. (a, b) the stress–strain curve, (c, d) the change in the degree of saturation with the mean net stress, (e, f) the change in the void ratio with the mean net stress.

35 / 46

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 9. The computed stress paths. (a) without hydraulic hysteresis, (b) with hydraulic hysteresis. 101

step #75 step #95 step #125

100

M

10-1 10-2 10-3

ED

rm rm  /rm  ext， n+1 int， n+1 2 ext， n+1 2

a

0

PT

10-5

102

5

10 15 20 25 30 35 40 45 50 55

m, number of global iteraions

c

10

0

10

AC

10-4 10-6 10-8

10

-10

10

-12

1

step #75 step #95 step #125 2

step #75 step #95 step #125

-4

10-5

i i Rp，  Rh，  n+1 1 n+1

10-2 i fn+1 (kPa2)

10-2 10-3

CE

b

10-4

10-6 10-7 10-8 10-9 10-10

3

4

10-11

i, number of local iteraions

1.0

1.5

2.0

2.5

3.0

i, number of local iteraions

Fig. 10. Convergence profiles of the proposed numerical method with hydraulic hysteresis for the isotropic loading and shearing test at the time steps of #75, #95 and #125. (a) the global residual force, (b) the residual of consistent condition, and (c) the residual for the hardening rule and flow 36 / 46

ACCEPTED MANUSCRIPT

flue.

a

b Load, PT=12 kPa PT

Load, PT

Water ifiltration, R

c

d

z1

Surface loading,PT 10 kPa

0

t

R

Water infiltration, R 14.4 mm/h

0

t

CR IP T

0.8

q=0

0.2

pw=0 0

t=0

AN US

Bottom surface t

 Δx Δy Δz 0.05m

z0

t>0

 Δx Δy Δz 0.025m  Δx Δy Δz 0.0125m

N: 189, E: 80 Coarse mesh

N: 6561, E: 5120 Fine mesh

M

N: 1025, E: 640 Medium mesh

Fig. 11. The computational model with coarse mesh (a, b), medium mesh (c) and fine mesh (d) for

ED

rain infiltration in a soil column: (a) t = 0, a surface load with 12 kPa is applied on the top surface; (b) t > 0, the surface load (PT = 12 kPa) immediately reduces to 10 kPa and then the water drains

PT

out from the bottom surface together with water infiltration into the top surface. N and E indicate the node number and element number (eight-node isoparametric brick element), respectively and

AC

CE

x, y, z are the mesh sizes in the x, y and z directions, respectively.

37 / 46

ACCEPTED MANUSCRIPT

0

0.8

a

4

b

z = 0.8 m 0.6

A

0.8

d

40

80 120 160 200 240 Time (min)

z = 0.2 m

0.6

-6

0.4

10-2

10-3

10-3

10-4

10-4

#75

10-5 10-6

0

h

10-1

10-2

20

30

40

10-6

80 120 160 200 240 Time (min)

0.0

0

40

80 120 160 200 240 Time (min)

Coarse mesh, modified Cws, Cas Medium mesh, modified Cws, Cas

#76

Fine mesh, modified Cws, Cas

0

20

40

60

80

m, number of global iteraions

ED

m, number of global iteraions

#76

Coarse mesh, original Cws, Cas

10-5 10

40

AN US

0

100

g

10-1

0.0

z = 0.2 m

#75

0.4

M

rm rm  /rm  ext， n+1 int， n+1 2 ext， n+1 2

100

80 120 160 200 240 Time (min)

80 120 160 200 240 Time (min)

0.8

0.2

40

40

f

1.2

-8 0

0

uz (mm)

pa (kPa)

-4

-12

1.6

e

z = 0.2 m

-2

pw (kPa)

-8

0

80 120 160 200 240 Time (min)

-4

CR IP T

0

40

0.2 0.0

A 0

uz (mm)

pa (kPa)

pw (kPa)

-12

0.4

z = 0.8 m

A

0

-4

-8

c

z = 0.8 m

Fig. 12. Evolutions of simulated results with different mesh sizes with respect to the center point at

PT

a height of (a, b, c) z  0.8 m and (d, e, f) z  0.2 m. (a, d) pore water pressure; (b, e) pore air pressure; (c, f) vertical displacement. (g)-(h) convergence profiles for the global residual force at

AC

CE

the time steps of #75 and #76. The time steps (#75 and #76) are labeled in (f).

38 / 46

ACCEPTED MANUSCRIPT

15

F

A

Rain intensity:16.0 mm/h

Z (m)

12 9

3

P

2

6

E

D

3

0

B

5

10

15

X (m)

20

CR IP T

0

5 25

30

8

35

C

AC

CE

PT

ED

M

AN US

Fig. 13. The finite element model for an homogenous soil slope under rain infiltration.

Fig. 14. The distributions of plastic deviatoric strain predicted by the coupled model without hydraulic hysteresis (a, c, e, g) and with hydraulic hysteresis (b, d, f, h) at t = 6 h (a, b), 7 h (c, d), 7.3 h (e, f), 7.5 h (g) and 8.4 h (h), respectively.

39 / 46

CR IP T

ACCEPTED MANUSCRIPT

Fig. 15. The distributions of the magnitude of displacement vector (a, b) and the degree of saturation (c, d) predicted by the coupled model without hydraulic hysteresis (a, c) and with

0.30

AN US

hydraulic hysteresis (b, d) at t = 7.5 h (a, c) and 8.4 h (b, d), respectively.

0.00

a

0.25

Uy (m)

without hysteresis with hysteresis

0.10 0.05 0.00

0.95

c

-0.20 0.7 0.5 0.4

pq (-)

0.3

AC

0.85 0.80

0

1

d

0.6

CE

0.90

-0.10 -0.15

PT

-0.05

Sr (-)

M

0.15

0.75

b

-0.05

ED

Ux (m)

0.20

0.2

A1 (6.88 h)

0.1 0.0 -0.1

2

3

4

5

6

7

8

9

-0.2

Time (h)

A2 (7.83 h) 0

1

2

3

4

5

6

7

8

9

Time (h)

Fig. 16. The evolutions of simulated results at the typical point P. (a) horizontal displacement; (b) vertical displacement; (c) the degree of saturation and (d) plastic deviatoric strain.

40 / 46

101

100

a Xmn+12/Xm+1 Xn2 n+1

rmext， n+1rmint， n+12/rmext， n+12

ACCEPTED MANUSCRIPT

100 10-1

t=3h t=6h t = 8.4 h (failure)

10-2 10

-3

0

20

40

60

80

100

10

b

-1

10-2 10-3 10-4 t=3h t=3h t = 8.4 h (failure)

10-5 10-6

20

40

60

80

m, number of global iteraions

CR IP T

m, number of global iteraions

0

Fig. 17. Convergence profiles of the proposed numerical method with hydraulic hysteresis at t = 3 h, 6 h and 8.4 h (slope failure). (a) the global residual force, and (b) the incremental value of the

AC

CE

PT

ED

M

AN US

primary variables.

41 / 46

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 18. The computed stress paths (a) and hydraulic paths (b) by the coupled model with/without

AC

CE

PT

ED

M

htdraulic hysteresis with respect to the point P.

42 / 46

ACCEPTED MANUSCRIPT

Tables (6 in total) Table 1. Integration procedure for the coupled hydro-mechanical model. 1. Initialize: p i0,  np,0 , pc (0)0n 1  pc (0)n , 0n1  0 , 1   n

2. Integration of hydraulic path: 2.1 Compute Se,trial n 1 , S e,w and S e,d

 

 

CR IP T

sn 1  sn  sn , en 1  en  en  n 1   n   n

 snkss1 mv exp mkss en1  v 1m S : 1  k m 1  S   e,n  snss v exp mkssv en   Se,w : [1  (  w exp( kp en1 ) sn1 )nv ] mv

 mv

AN US

trial e,n 1

Se,d : [1  (  d exp( kp en1 ) sn1 )nv ] mv

M

2.2 Check

trial If Se,w  Se,trial n 1  Se ,d , then Se,n 1  Se,n 1

ED

Else if Se,trinal1  Se,d , then Se,n 1  Se,d

PT

Else Se,n 1  Se,d

CE

3. Elastic predictors:

 ni1 : Ce1 :   n 1   np,i1 

AC

4. Check the residuals:

 n1 

1  S r,1/4n1 1 i ; pni1    iii,n1 ; qni 1   3 J 2 n1 ; g ( en ) 3

 (   )ln pc (0)in1  N ( h( n1 )  1)  pc ( )in1  exp   h( n1 )    

43 / 46

ACCEPTED MANUSCRIPT f ni1 :  qni 1   M 2  pni1   pc ( )in 1  pni1  2

2

i

i p,n 1

: 

i h,n 1

:  pc (0)

R

R

p,i n 1

    p n

i n 1

i n 1

g   n 1

i  1 e g  i  pc (0)n exp   n 1      n 1   p  

5. Compute:

CR IP T

i i i If f n 1   f and Rp,n 1 1  Rh,n 1   ph then Exit; else Goto 5

f g , g , f ,   n 1 p  n 1   n 1 pc (0) i

i

i

i

i

i

n 1

i

2 g 2 g 2 g , , ,  pc (0) n 1 p  n 1 ppc (0) n 1 n 1

AN US

2 g  2

i

i

 1 e g 1  e    pc (0)n exp   p     n 1   

ED

M

i i   2 g 2 g 1 i i  C e   n1  n 1 2         p (0)   c i n 1 n 1 1 i  n1 :   ,  i i n 1 2 g 2 g   i i  n1  1   n1  p  ppc (0) n1  n 1 

PT

6. Obtain the increment of plastic multiplier:

AC

CE

 2  in 1 

f

i n 1

 f    

 f     f  pc (0) 

f  pc (0) 

i

 Ri   1  p,i n +1   Rh,n +1  n 1

i

 g g   1   p     n 1

T i

n 1

7. Obtain the incremental plastic strain and yield stress:

 g i      n 1   Ri  2 i   ni1  1  p,n +1  1      i     n 1   i  i    p (0)  Rh,n +1   c n 1   g   p   n 1 

8. Update state variables and plastic multiplier:

44 / 46

ACCEPTED MANUSCRIPT  ni11   ni1   ni1  np,i11   np,i1  C e1 :  ni1  in11   in 1   2  in 1 pc (0)in11  pc (0)in 1  pc (0)in 1

CR IP T

9. i i 1, Goto 3.

Table 2. Model parameters for the coupled hydro-mechanical constitutive model. Mechanical model Value

Parameter

N

1.52 0.13 0.02 0.72 3000 MPa 30.2 kPa 6.1427 0.8123

nw mw kp

d w

kss kse Sr, res

ED

M G pc(0) a b

M

 

AN US

Parameter

Water retention curve mode Value Without hysteresis With hysteresis 0.7992 0.7992 0.1554 0.1554 10.605 10.605 3.1120107 kPa-1 1.3469 106 kPa-1 1.7038106 kPa-1 0.072 0.857 0 0

Table 3. Material properties for the soil column.

AC

CE

PT

Parameters Bulk modulus of solid grains Ks (MPa) Solid grain density s (kg/m3) Initial porosity 0 Initial intrinsic permeability k0 (m2) Bulk modulus of water Kw (MPa) Initial water densityw0 (kg/m3) Water viscosity w (kg/ms) Initial air densitya0 (kg/m3) Air viscositya (kg/ms) Atmospheric ref. pressure patm (kPa)

45 / 46

Value 1.0106 2850 0.2975 4.51013 2.0103 1000 1.0103 1.25 1.8105 101.0

ACCEPTED MANUSCRIPT

Table 4. Model parameters for the coupled hydromechanical constitutive model.

  M G a b -

Water retention curve mode Parameters Value nw 6.6787 mw 0.2525 kp 1.6105 d 0.0515 kPa1 w 0.0800 kPa1 kss 0.38 kse 0.58 Sr, res 0

Value 1.62 0.2 0.05 1.0 3000 MPa 10.0 1.5 -

CR IP T

Mechanical model Parameters N

Table 5. Comparisons of the simulated results at the end of simulation with the coarse, medium

–2.0598 0.0132 0.8317 –6.0029 0.1056 –8.0110

–2.0609 0.0130 0.8456 –5.9804 0.1027 –7.9909

PT

z = 0.8 m

pw (kPa) Pa (kPa) uz (mm) pw (kPa) Pa (kPa) uz (mm)

Coarse mesh Medium mesh (Vc) (Vm)

M

z = 0.2 m

Variables

ED

Height

AN US

and fine meshes.

Fine mesh (Vf) –2.0613 0.0129 0.8483 –5.9778 0.1025 –7.9802

Relative differences (VcVf)/Vf –0.07% 2.33% –1.96% 0.42% 3.02% 0.39%

(VmVf)/Vf –0.02% 0.78% –0.32% 0.04% 0.20% 0.13%

AC

CE

Table 6. Model parameters for the coupled hydromechanical constitutive model. Mechanical model Water retention curve mode Value Parameter Value Parameter Without hysteresis With hysteresis N 1.32 nw 1.56 1.560 0.143 mw 0.36 0.361  0.023 kp 2.21 2.212  M 1.0 d 0.02857 kPa1 0.02857 kPa-1 G 3000 MPa w 0.07931 kPa1 a 9.563 kss 0.38 b 1.910 kse 0.58 Sr, res 0 0

46 / 46