# Analysis of cement-treated clay behavior by micromechanical approach

Dong-Mei Zhang, Zhenyu Yin, Pierre-Yves Hicher, Hong-Wei Huang

► To cite this version:

Dong-Mei Zhang, Zhenyu Yin, Pierre-Yves Hicher, Hong-Wei Huang. Analysis of cement-treated clay behavior by micromechanical approach. Frontiers of Structural and Civil Engineering, 2013, 7 (2), pp.137-153. 10.1007/s11709-013-0204-z. hal-01005247

HAL Id: hal-01005247

https://hal.science/hal-01005247

Submitted on 18 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Analysis of cement-treated clay behavior by micromechanical approach

Dong-Mei ZHANGª, Zhen-Yu YINᵇ, Pierre-Yves HICHERª,ᶜ, Hong-Wei HUANGª

**ABSTRACT** Experimental results show the significant influence of cement content on the mechanical properties of cement-treated clays. Cementation is produced by mixing a certain amount of cement with the saturated clay. The purpose of this paper is to model the cementation effect on the mechanical behavior of cement-treated clay. A micromechanical stress-strain model is developed considering explicitly the cementation at inter-cluster contacts. The inter-cluster bonding and debonding during mechanical loading are introduced in two ways: an additional cohesion in the shear sliding and a higher yield stress in normal compression. The model is used to simulate isotropic compression and undrained triaxial tests under various confining stresses on cement-treated Ariake clay and Singapore clay with various cement contents. The applicability of the present model is evaluated through comparisons between numerical and experimental results. The evolution of local stresses and local strains in inter-cluster planes are discussed in order to explain the induced anisotropy due to debonding at contact level under the applied loads. The numerical simulations demonstrate that the proposed micromechanical approach is well adapted for taking into account the main physical properties of cement-treated clay, including damage and induced anisotropy under mechanical loading.

**KEYWORDS** microstructure, cementation, clay, micromechanics, anisotropy, debonding

## 1 Introduction

Stabilization of soft ground by deep cement mixing and jet grouting methods, which both introduce cement into the ground, has been widely used in geotechnical projects [1–4]. Due to the increase of their mechanical properties, the cement-treated soil layer helps to decrease the soil movements during construction. Thus, proper understanding of the cementation influence on the mechanical behavior of cement-treated clay is of great importance.

The improvement of the engineering properties of cement-treated soil has been attributed to the soil-cement reaction [5,6]: Hydration reaction and pozzolanic reaction between the hydrated lime and the clay minerals. This creates a type of microstructure different from the natural

structured soils (soft sensitive clay and stiff clay [7,8]), for which soil structures are often formed by solute deposition at inter-particle contacts, charge deficiencies, and Van der Waal's forces [9].

Experimental results show an increase of shear strength and shear modulus with cement content, accompanied by a large post-peak reduction in the strength of the treated soil [3–5,10,11]. Cementation also makes the yield stress and swelling index higher than those of untreated soil. Bond degradation is often observed if the compression load exceeds the yield stress, as it has also been observed in natural structured soils [12–14]. The difference between natural structured clay and cement-treated clay lies in the properties of host clay modified by adding cement, which needs to be considered in constitutive modeling.

The behavior of cement-treated clay has been modeled by conventional elasto-plastic models, e.g., Lee et al. [15],

<footnote>ªDepartment of Geotechnical Engineering, Tongji University, Shanghai 200092, China</footnote>

<footnote>ᵇDepartment of Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China</footnote>

<footnote>ᶜResearch Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale de Nantes, BP 92101, Nantes 44321, France</footnote>
Horpibulsuk et al. [16]. In their approaches, the material structure, including fabric and bonding, was accounted for by the initial size of the yield surface. A damage-type mechanism was introduced which produced the reduction of the yield surface size due to bond degradation. However, the bonding and debonding at inter-cluster contact level has not been explicitly considered in these models.

To tackle this problem, a micromechanical modeling appears to be a suitable tool. The approach is based on a homogenization technique in order to derive the stress-strain relationship of the granular assembly from forces and displacements at the particle level [17–29]. Along these lines, Chang and Hicher [27] have proposed a model for granular materials based on the development of micro–macro links between strain and inter-particle displacements on one hand, and stress and inter-particle forces on the other hand, which allow both the micro-system (inter-granular forces and displacements) and the macro-system (overall stress–strain) to be in balance. More recently, the model has then been extended to clays by considering the clayey microstructure as a collection of clusters which play the same role as the grains in granular assemblies for various features [7,8,30–33].

Horpibulsuk et al. [3] have shown that the clay fabric consists of clay cluster (aggregation of clay particles) and that the role of the induced cementation is to weld clay clusters (see Fig. 1(a)). Kamruzzaman et al. [4] have made SEM photos of untreated and cement treated Singapore clay with various cement contents under various con-solidation pressures (see Fig. 1(b)), which demonstrate that the size of the clay clusters (as shown by dash circle in Fig. 1(b)) can be considered the same for all cases. Therefore, the cement-treated clay can be seen as an assembly of clusters (aggregates of clay-cement particles). The deformation of an assembly can be obtained by integrating the movement of the inter-cluster contacts in all orientations. Thus, the effect of the inter-cluster bonding can be explicitly taken into account. As the orientation dependent properties are explicitly introduced, the induced anisotropy due to inter-cluster debonding can be modeled

in a direct way.

In this study, the cementation effect on the mechanical properties of cement-treated clays is first investigated. Based on this, the dependency of the inter-cluster bonding and of the critical state line (CSL) location on the cement-content are implemented into the model developed by Yin et al. [7] for sensitive clay. The extended model is then used to simulate isotropic compression and undrained triaxial tests under various confining stresses on cement-treated clays with different cement contents. The evolution of local stresses and strains in inter-cluster planes due to externally applied loads is also discussed in order to describe the mechanism of induced anisotropy due to debonding at the micro contact level. The overall applicability of the present model is evaluated based on the comparison between experimental results and numerical simulations.

## 2 Influence of cement content on mechanical properties

In this section, the cement-treated Ariake clay with cement content varying from 0% to 18% (0% representing untreated clay [3]) and the cement-treated Singapore clay with cement content varying from 0% to 50% [4] were selected to investigate the influence of cement content on clay mechanical properties. The percentage of cement content was calculated based on the ratio of dry weight of cement to dry weight of soil. Some physical properties of the untreated clays are summarized in Table 1. The type I ordinary Portland cement was used for treating both clays.

The initial water content of the clay slurry may influence the initial void ratio, the compression yield stress, and the destructuring rate during mechanical loading for samples of cement-treated clay. In this paper, the two clay slurry were prepared at a water content equal to 1.4–1.5 times the liquid limit (180% for Ariake clay, 120% for Singapore clay), and the influence of the initial water content was not studied. Furthermore, all the experimental results selected

<figure><img src="image_2.png" /><figcaption>**Fig. 1** (a) Microfabric of uncemented clay and structure of induced cemented clay (after Horpibulsuk et al. [3]); (b) SEM images of 10% and 50% cement-treated clay at 28-day curing (after Kamruzzaman et al. [4])</figcaption></figure>
<table><caption>**Table 1** Index properties of Ariake clay and Singapore clay</caption><thead><tr><th>index property</th><th>Ariake clay</th><th>Singapore clay</th></tr></thead><tbody><tr><td>water content, w/%</td><td>135–150</td><td>70</td></tr><tr><td>liquid limit, wL/%</td><td>120</td><td>87</td></tr><tr><td>plastic limit, wP/%</td><td>57</td><td>35</td></tr><tr><td>clay particle size content/%</td><td>55</td><td>68</td></tr><tr><td>silt content/%</td><td>44</td><td>22</td></tr><tr><td>sand content/%</td><td>1</td><td>10</td></tr><tr><td>specific gravity/G<sub>s</sub></td><td>2.60</td><td>2.67</td></tr></tbody></table>

for this study were obtained on samples of cement-treated clay with a curing time of 28 days. The influence of curing time on the cementation was not studied at current stage.

## 2.1 Cement-content dependent behavior under isotropic and 1D compression

### Isotropic compression and oedometer tests by Horpibulsuk

et al. [3] and Kamruzzaman et al. [4] on two cement-treated clays show the significant influence of the cement content on the swelling index (see Fig. 2(a)) and on the compression yield stress (see Fig. 2(b)). This influence can be expressed by Eq. (1) for the swelling index $κ$ and by Eq. (2) for the bonding pressure $p_b$ defined as the difference of the compression yield stress between cement-treated clay $p_{yc}$ and untreated clay $p_{yu}$ : $p_{yb} = p_{yc} - p_{yu}$, as follows:

$$ \kappa = \kappa_0 \exp(-\beta_k c), \quad (1) $$

$$ p_b = \beta_{p1} [\exp(\beta_{p2} c) - 1], \quad (2) $$

where $κ_0$ is the swelling index for untreated clay; c is the cement content; $β_k$, $β_{p1}$ and $β_{p2}$ are material constants (see Fig. 2(a), (b)).

By comparing the post-yield stress-strain curves of cement-treated samples to those of untreated samples [3,4], a mechanical bond degradation process can be observed. This bond degradation is an irreversible phenomenon that,

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Effect of cement content on properties of cement treated clays. (a) Swelling index; (b) bonding pressure; (c) friction angle at failure, and (d) apparent bonding cohesion</figcaption></figure>
experimentally, appears to be controlled by plastic strain accumulation [34]. The bonding and debonding can be linked directly to the inter-cluster bonds due to the cementation formed when adding cement into clay slurry.

## 2.2 Cement-content dependent shear behavior

The friction angle at failure was measured based on undrained triaxial tests by Horpibulsuk et al. [3] and Kamruzzaman et al. [4] on both cement-treated clays by extending the curves to an axial strain up to 15% according to the previous trend. As shown in Fig. 2(c), for Ariake clay, having an initially high friction angle, the cement content influences slightly the friction angle at failure (we obtained $\phi_{\mu} = 40^{\circ}$ measured from undrained stress paths, which is slightly different from $\phi_{\mu} = 38^{\circ}$ provided by Horpibulsuk et al. [3]). However, for Singapore clay with an initially small friction angle, the cement content influences significantly the friction angle at failure. This influence can be expressed by Eq. (3) as follows:

$$ \phi_{\mu} = \phi_{\mu 0} + (\phi_{\mu max} - \phi_{\mu 0})[1 - \exp(-\beta_{f}c)], \quad (3) $$

where $\phi_{\mu 0}$ is the friction angle at failure for the untreated clay; $\phi_{\mu max}$ is the maximum friction angle at failure for the cement-treated clay; $\beta_{f}$ is a material constant.

The apparent cohesion due to cement bonds is measured from the position of the failure line in the $p'-q$ plane. Figure 2(d) shows the significant influence of the cement content on the apparent cohesion $c_b$ for both cement-treated clays. This influence can be expressed by Eq. (4), having a form similar to Eq. (2) for the bonding pressure, as follows:

$$ c_b = \beta_{c1}[\exp(\beta_{c2}c) - 1], \quad (4) $$

where $\beta_{c1}$ and $\beta_{c2}$ are material constants.

The stress-strain curves of the two cement-treated clays show also the influence of the cement content on the initial shear modulus, which is similar to its influence on the swelling index (or bulk modulus). Since the shear modulus cannot be accurately measured based on digitized data, we will not discuss it in this paper.

Upon increase of the applied shear stresses, the inter-cluster bonds are progressively damaged. As a result, a decrease of the deviatoric stress takes place after the peak for cement-treated samples.

## 2.3 Cement-content dependent critical state line (CSL)

The critical state lines are obtained from the undrained triaxial tests for untreated clays, as shown in Figs. 3(a) and (b). The void ratio corresponding to this state is $e_c$. The CSL can be written as follows for clay:

$$ e_c = e_{cr} - \lambda \ln \left( \frac{p'}{p_{ref}} \right), \quad (5) $$

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3 Effect of cement content on the location of critical state line. (a) For cement treated Ariake clay; (b) for cement treated Singapore clay; (c) reference critical void ratio versus cement content</figcaption></figure>
where $p'$ is the mean effective stress; $\lambda$ is the slope of the critical state line in the e-logp' plane, $e_{cr}$ is the reference void ratio at critical state corresponding to $p_{ref} = 100$ kPa.

For cement-treated clays, due to the strain localization of samples during shearing, the critical state line is difficult to be obtained accurately. If the failure state in the $p'^{-}q$ plane (corresponding to 15% of axial strain) is considered as the critical state, the critical state lines in the e-logp' plane can be obtained for all cement-treated samples. For cement-treated Ariake clay (see Fig. 3(a)), the critical state line can be assumed to be parallel to that of the untreated clay. If this assumption is also adopted for cement-treated Singapore clay, the critical state lines can be obtained as shown in Fig. 3(b). Therefore, the influence of cement content on the location of the critical state line can be obtained by making the reference void ratio function of the cement content, as shown in Fig. 3(c), expressed as follows

$$e_{cr} = e_{cr0} + (e_{crmax} - e_{cr0})[1 - \exp(-\beta_{c}c)], \quad (6)$$

where $e_{cr0}$ is the reference critical state void ratio for the untreated clay; $e_{crmax}$ is the maximum reference critical state void ratio for cement-treated clay; $\beta_{c}$ is a material constant.

# 3 Micromechanical model

The model developed in this study is an extension of the model developed by Yin et al. [7] for sensitive clay. The modifications made to account for the effect of cement content on the inter-cluster bonding and debonding and on the CSL location are discussed below. As stated in the previous section, the bonding and debonding will be considered for two different loading conditions: normal compression and shear sliding.

## 3.1 Inter-cluster behavior

To take into account the effect of the inter-cluster bonding, the plastic law of the shear sliding at inter-cluster contacts proposed by Yin et al. [7] is adopted. The yield function for shear sliding is assumed to be of Mohr-Coulomb type, written in a contact plane (e.g., $\sigma$, $\tau_{s}$, $\tau_{t}$ defined in Appendix A, see coordinate system in Fig. 4) as follows

$$F_1(\sigma, \tau, c_b, H_1) = \frac{\tau}{\sigma + c_b / \tan \phi_\mu} - H_1(\gamma^p), \quad (7)$$

where $c_b$ is the inter-cluster cohesion due to the inter-cluster bonding, $H_1(\gamma^p)$ is a hardening/softening function. The values of $c_b$ and $\phi_\mu$ depend upon the cement content (see Eqs. (3) and (4)). The influence of the cement content on the shear behavior is thus introduced. For untreated clay, $c_b = 0$, and the equation can be reduced to that used by Yin et al. [33].

The local dilatancy equation is derived from that

proposed by Yin et al. [33], as follows

$$\frac{d\varepsilon^p}{d\gamma^p} = D \left| \left( \tan \phi_\mu - \tan \phi_m \right) \left( \frac{\tan \phi_m}{\tan \phi_\mu} \right) \right| \left( \frac{e}{e_c} - 1 \right), \quad (8)$$

where $\tan \phi_m = \tau / (\sigma + c_b / \tan \phi_\mu)$; D is a material constant which controls the dilatancy rate during shearing. The void ratio at critical state $e_c$ is calculated using Eq. (5) where the influence of the cement content is introduced by $e_{cr}$ (see Eq. (6)).

The degradation of the inter-cluster bonds can be modeled as a damage of the bonded contacts. Therefore, the damage law proposed by Yin et al. [7] is adopted in the expression of the inter-cluster cohesion, as follows

$$c_b = c_{b0} \exp(-\xi_d \gamma^p), \quad (9)$$

where $c_{b0}$ is the initial inter-cluster cohesion; $\xi_d$ is the factor of damage representing the influence of the tangential displacement in the damage law. Therefore, during loading each contact sliding produces a progressive damage of the cohesion.

To account for the inter-cluster bonding effect on the compressive behavior of cement-treated clay, we adopt the plastic law for the normal compression of inter-cluster contacts proposed by Yin et al. [7]

$$F_2(\sigma, \sigma_p) = \sigma - \sigma_p, \quad (10)$$

where $\sigma_p = \sigma_{pi}(1 + \chi)$ is the compression yield stress; $\sigma_{pi}$ is the intrinsic compression yield stress corresponding to the yield stress measured from an isotropic compression test on a reconstituted sample of untreated clay; $\chi$ is the bonding ratio defined by $\chi = \sigma_p / \sigma_{pi} - 1$. The value of $\sigma_p$ depends on the cement content ($\sigma_p = \sigma_{pi} + p_b$, see Eq. (2) for $p_b$). Thus, the influence of the cement content on the compression behavior is introduced.

The hardening function controlling the evolution of $\sigma_{pi}$ is defined as follows

$$\sigma_{pi} = \sigma_{pi0} \exp \left( \frac{\varepsilon^p}{c_p} \right), \quad (11)$$

where $c_p$ is the compression index determined from the compression curve plotted in the $\varepsilon^p - \log \sigma$ plane.

The damage law proposed by Yin et al. [7] is also adopted in the expression of the bonding ratio at inter-cluster level,

$$\chi = \chi_0 \exp(-\xi_n \varepsilon^p), \quad (12)$$

where $\chi_0$ is the initial bonding ratio depending on the cement content; $\xi_n$ is the factor of damage representing the influence of the normal displacement in the damage law. Therefore, during loading each contact produces a progressive damage of the bonding ratio. For untreated clay $\chi_0 = 0$ (or $\sigma_p = \sigma_{pi}$), the equation can be reduced to that used by Yin et al. [33].
<figure><img src="image_7.png" /><figcaption>Fig. 4 Local coordinate at inter-particle contact</figcaption></figure>

## 3.2 Study of cementation related parameters

The parametric study of the cementation related parameters of the model is carried out under different loading conditions using the model constants given in Table 2 for the cement-treated Ariake clay. First, an isotropic compression test with $p_b = 25$ kPa is considered. The numerical simulations using different values of the damage parameter $\xi_n$ are shown in Fig. 5(a). In this figure, the destructuration process is controlled by the destructuration rate which depends on the damage parameter $\xi_n$. Based on this, different values of $\xi_n$ can be determined for different amounts of cement.

Then, an undrained triaxial test on a normally consolidated sample is selected for the simulation, as shown in Figs. 6(a) and (b) for $\xi_d = 1$ and different $c_{b0}$ (20, 65 and 140 kPa), and in Figs. 6(c) and (d) for $c_{b0} = 140$ kPa and different $\xi_d$ (from 0 to 10). The initial bonding cohesion $c_{b0}$ controls the magnitude of the peak strength whereas $\xi_d$ controls the destructuration rate of the cohesion. Figures 6(e) and (f) show the coupling effect

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 5 Compression behavior. (a) Effect of $\xi_n$ on destructuration; (b) cement-treated Ariake clay; (c) cement-treated Singapore clay</figcaption></figure>

of $c_{b0}$ and $\xi_d$. Different values of ($c_{b0}$, $\xi_d$) can be selected in order to obtain the same stress level at axial strain equal to 15%, but with different stress evolutions during straining.
<table><caption>**Table 2** Values of model parameters for cement-treated clays</caption><thead><tr><th>parameters</th><th>e<sub>0</sub></th><th>κ<sub>0</sub></th><th>λ</th><th>σ<sub>pit0</sub></th><th>φ<sub>μ0</sub></th><th>k<sub>pR</sub></th><th>D</th><th>e<sub>cr0</sub></th><th>ξ<sub>n</sub></th><th>ξ<sub>d</sub></th></tr></thead><tbody><tr><td>Ariake clay</td><td>3.44</td><td>0.086</td><td>0.446</td><td>10</td><td>40°</td><td>0.15</td><td>6</td><td>1.911</td><td>10</td><td>3</td></tr><tr><td>Singapore clay</td><td>2.50</td><td>0.092</td><td>0.355</td><td>10</td><td>23°</td><td>1.00</td><td>10</td><td>1.500</td><td>20</td><td>5</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 6** Parametric study for undrained shearing behavior. (a) and (b) effect of $c_b$; (c) and (d) effect of $ξ_d$; (e) and (f) coupling effect for $c_b$ and $ξ_d$</figcaption></figure>

Therefore, $c_{b0}$ and $ξ_d$ can be obtained by fitting the experimental stress-strain curves of the cement-treated clay.

# 4 Test simulations

In this section, the model is used to simulate isotropic compression and undrained triaxial tests under various confining stresses on untreated and cement-treated Ariake clay and Singapore clay with different cement contents. All undrained triaxial tests were simulated by taking into account the first stage of isotropic compression before shearing as imposed on the specimen in the laboratory. The applicability of the model is investigated by comparing its predictions with experimental results.

## 4.1 Determination of model parameters

### 4.1.1 Ariake clay

Based on the isotropic compression and undrained triaxial tests under a consolidated stress of 200 kPa on the untreated Ariake clay (see Fig. 5(b) “c = 0%”, Fig. 7-1 “cement content: 0%”), several parameters can be determined (see Table 2).

1) Contact number per unit volume $N/V$ and mean aggregate size $d$. Based on SEM photos, the mean size $d$ of the clay aggregates is assumed to be 10 µm, and $N/V$ is obtained by the mean size $d$ and the void ratio $e$.

2) Inter-cluster elastic constants: $k_{n0}$, $k_{pR}$ and $n$. The exponent $n = 1$ is usually assumed in order to obtain a linear rebound line in the $e-\log p'$ plane. $K_{pR} = 0.5$ can
usually be assumed for clay (see [7,8,30–33]). The contact
stiffness $k_{n0}$ can be calculated by using the swelling index
$\kappa$, measured from the swelling curve by using the relation
$k_{n0} = 9(1 + e_0)\sigma_{at}/(4\kappa)$. Therefore, only $\kappa_0$ and $e_0$ are
needed as input parameters. For the Ariake clay, $\kappa_0 = 0.086$
and $e_0 = 3.44$ were measured.

3) Inter-cluster normal compression: $\sigma_{pi0}$ and $c_p$. The initial value of the intrinsic compression yield stress $\sigma_{pi0} = 10$ kPa is assumed for the untreated Ariake clay, since there is no measurement available. Using the intrinsic compression index $\lambda_i$ and the swelling index $\kappa$ of an assembly of clay aggregates, we can also determine the plastic inter-cluster compression index $c_p$ by the relation $c_p = (\lambda_i - \kappa)/3(1 + e_0)$. Therefore, $\lambda_i$ is used as an input parameter instead of $c_p$. For the Ariake clay, $\lambda_i = 0.446$ was measured.

4) Inter-cluster friction angle: $\phi_{\mu}$ and $m$. $m = 1$ is generally assumed for clay (see [7,8,30–33]). $\phi_{\mu0}$ is assumed equal to the macroscopic friction angle of the normally consolidated reconstituted material. For Ariake clay, $\phi_{\mu0} = 40^{\circ}$ was measured.

5) Inter-cluster plastic shear stiffness ratio $k_{pR}$. According to the hardening rule used for shear sliding, a bigger value of $k_{pR}$ provides a stiffer shear stress-strain curve. For Ariake clay, $k_{pR} = 0.15$ was obtained by fitting the deviatoric stress versus axial strain curve.

6) Dilatancy constant $D$. According to the local dilatancy relation, a bigger value of $D$ gives bigger plastic volumetric changes, which results in higher excess pore pressure in undrained condition. Based on this, $D = 6$ was obtained for Ariake clay by fitting the effective stress path of the selected undrained test.

7) Critical state for the soil $e_{cr0}$. $e_{cr0}$ is the reference void ratio at critical state corresponding to $p_{ref} = 100$ kPa. For Ariake clay, $e_{cr0} = 1.911$ was measured.

8) Cementation related parameters. Based on the analysis of the influence of cement content on the mechanical properties of cement-treated Ariake clay (see Section 2, Figs. 2 and 3), the cementation related parameters are summarized.

9) Inter-cluster debonding related parameters $\xi_n$ and $\xi_d$. Based on two selected isotropic compression and undrained triaxial tests under a consolidated stress of 200 kPa on cement-treated Ariake clay (see Fig. 5b "c = 6%", Fig. 7(2) "cement content: 6%"), debonding related parameters can be determined (see Table 2). $\xi_n$ controls the debonding rate of $\chi$. For a given compression curve $\xi_n = 10$ was obtained (see discussion in Section 3.2). The parameter $\xi_d = 3$ was obtained by curve fitting of the selected undrained triaxial test (see discussion in Section 3.2).

All determined values of the parameters with the cement-dependent parameters (Figs. 2 and 3) were used to predict the other tests.

### 4.1.2 Singapore clay

The determination of parameters for cement-treated Singapore clay were based on the analysis of the influence of cement content on the mechanical properties of cement-treated Singapore clay (see the section 2, Figs. 2 and 3), considering one one-dimensional test (see Fig. 5(c) "c = 0%"), one isotropic compression test (see Fig. 5(c) "c = 10%"), and one undrained triaxial test under a consolidated stress of 300 kPa on cement-treated Singapore clay (see Fig. 8-1 "cement content: 10%"). The procedure of parameter determination is similar to the one described for Ariake clay and is not repeated here. All determined parameters, summarized in Table 2, were used to predict all the other tests.

## 4.2 Numerical simulations of isotropic and 1D consolidation

Figure 5(b) shows the comparison between experimental results and numerical simulations of isotropic compression tests on untreated and cement-treated Ariake clay with different cement contents. Figure 5(c) shows the same comparison for one-dimensional compression test on untreated Singapore clay and isotropic compression tests on cement-treated Singapore clay with different cement contents. For both clays, a satisfactory agreement is achieved by using for each clay the corresponding set of parameters in Table 2.

## 4.3 Numerical simulations of undrained triaxial tests

Figure 7 shows the comparisons for all the selected undrained triaxial tests on untreated and cement-treated Ariake clay with different cement contents by using for all the numerical simulations the same set of parameters presented in Table 2. Similarly, Fig. 8 shows the comparisons for all selected undrained triaxial tests on cement-treated Singapore clay with different cement contents. In the whole, the numerical simulations agree with the experimental results, although discrepancies were found for some of the tests, which may be due to some variations of the sample physical characteristics. The model appears to be able to reproduce the mechanical behavior of cement-treated clay and the effects due to the cementation and the damage process during shearing with reasonable accuracy.

Figure 9 shows the comparisons between experiments and simulations for undrained triaxial tests on samples of Ariake clay with different cement contents consolidated up to 200 kPa. The influence of the cement content on the behavior of cement-treated Ariake clay was well reproduced by the proposed model: 1) The increase of stiffness and shear strength with cement content, accompanied by a
<figure><img src="image_12.png" /><figcaption>**Fig. 7** Comparisons between experiments and simulations for undrained shearing tests on cement treated Singapore clay with cement content of (1) 0%; (2) 6%; (3) 9%; (4) 12%; (5) 18%. (a) Effective stress path; (b) deviatoric stress versus axial strain; (c) excess pore pressure versus axial strain</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 8** Comparisons between experiments and simulations for undrained shearing tests on cement treated Singapore clay with cement content of (1) 10%; (2) 30%; (3) 50%. (a) Effective stress path; (b) deviatoric stress versus axial strain; (c) excess pore pressure versus axial strain</figcaption></figure>

large post-peak reduction in the strength of the treated soil; 2) samples with low cement content (c = 0%, 6%, 9%, 12%) exhibit a contractive behavior and samples with high cement content (c = 18%) a dilative one. Figure 10 shows the comparisons between experiments and simulations for undrained tests on Singapore clay with different cement contents consolidated up to 500 kPa. Again, the influence of cement content was well reproduced by the proposed model.

# 5 Micromechanical analysis of induced anisotropy

In this section, we investigate the local stress-strain

behavior for selected contact planes with different orientations θ = 0°, 18°, 28°, 45°, 55°, 72° and 90° (see Fig. 4, x representing vertical direction). The undrained triaxial tests on samples of Ariake clay with cement contents of 0%, 9% and 18% under a consolidation stress of 200 kPa (see Fig. 9) were chosen as an example for illustrating the local behavior in contact planes (see Figs. 11 and 12). The influence of the cement content on the local stress-strain behavior is thus discussed. To study the stress-induced anisotropy during undrained shearing, we selected for all the tests three steps for plotting the variable distribution in the rose diagrams (see Fig. 13): the beginning of loading, the peak deviatoric stress and the end of loading.
<figure><img src="image_14.png" /><figcaption>**Fig. 9** Effect of cement content on undrained shearing behavior of cement treated Ariake clay for tests under confining pressure of 200 kPa. (a)–(c) Experimental results, (d)–(f) calculated results</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 10** Effect of cement content on undrained shearing behavior of cement treated Singapore clay for tests under confining pressure of 500 kPa: comparisons of (a) effective stress path; (b) deviatoric stress versus axial strain; (c) excess pore pressure versus axial strain</figcaption></figure>

## 5.1 Stress-strain behavior and bond degradation at local scale

Figure 11(a) presents the local stress paths for all selected contact orientations (from 11-1 to 11-5) for all cases (cement contents of 0%, 9% and 18%). For the 0° and 90° contact planes, the shear stress is null. Thus, the stress-strain behaviors are not plotted. For the others, the stress state with the highest stress ratio is located in the 70° contact plane. Figure 11(b) shows the local shear stress-strain curves which clearly indicate that every contact plane is mobilized to a different degree. The planes with

the largest movements are located around the orientation of 72° (close to $π/4 + φ_μ/2 = 65°$). These active contact planes contribute largely to the overall deformation of the specimen. Among all selected contact planes, only the 72° contact plane has a local shear stress-strain curve similar to the global one. Figure 11(c) presents the local normal stress-strain curves, in which normal strains are mostly composed by elastic strains due to normal stresses and plastic strains induced by the dilatancy law during shearing. The results show that, although the global volumetric strain is null during undrained shearing, local normal strains develop in every contact plane to different
<figure><img src="image_16.png" /><figcaption>**Fig. 11** Local stress-strain behavior for selected tests. (a) Stress path in shear stress versus normal stress; (b) shear stress versus shear strain; (c) normal strain versus normal stress</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 12** Local bond degradation for selected test. (a)-(b) shear bond cohesion versus global axial strain, and (c) normal bonding ratio versus global axial strain</figcaption></figure>

degrees and with different signs (negative or positive value). For each sample with different cement content, each contact plane behaves differently, and the influence of the cement content on the local stress-strain behavior is significant, which agrees with the global behavior (see Fig. 9).

The degradation of the bonding ratio and of the bonding cohesion for all selected contact planes are plotted versus the global axial strain in Fig. 12 for samples with cement contents of 9 and 18%. For the bonding cohesion (Figs. 12 (a)–(b), most of the degradation took place in the 72° contact plane due to its largest local displacements. The other contact planes have only a slight debonding because their mobilized displacements are small. For the bonding ratio (Fig. 12(c)), a very small reduction took place because the maximum normal strain remained smaller than 1% (see Figs. 11(c)).

## 5.2 Local stresses, strains and bond distributions

The distributions of the local normal stiffness and bond in contact planes of various orientations are plotted in Fig. 13 (Fig. 13-1 for untreated clay, Fig. 13-2 for cement content of 9% and Fig. 13-3 for cement content of 18%).

1) The normal stiffness is a function of the normal stress and of the amount of cement. The normal stiffness distribution at the end of the isotropic consolidation has a circular shape (see the bold line in Figs. 13(1-a), (2-a) and (3-a) which corresponds to an isotropic distribution of the normal stress for all contact planes. During the undrained shearing, different cement contents give different evolutions of the stiffness distribution: for 0% cement content, the distribution becomes anisotropic with an elliptical shape having the long axis in the vertical direction; for 9% cement content, the distribution becomes also anisotropic with higher stiffness in the vertical direction and lower stiffness in the horizontal direction, and then shrinks

during shearing due to bond damage; the same evolution is found for 18% cement content with higher stiffness values due to stronger bonds.

2) During the undrained shearing, the shear stress distribution expands while maintaining a similar shape for all cases (Figs. 13(1-b), (2-b) and (3-b)). This distribution shows a concentration of the maximum local shear stresses between the 45° and 55° orientations, in agreement with the local stress-strain curves in Fig. 11.

3) During the undrained shearing, the shear strain distribution for all cases (Figs. 13(1-c), (2-c) and (3-c)) shows that very large strains occurred at the end of the loading in contact planes around the 70° orientation, which agrees with the local stress-strain curves in Fig. 11.

4) The distribution of shear bonding cohesion changes slightly from step 1 to 2 due to small mobilized shear plastic strains for both 9% and 18% cement contents. Then, from step 2 to 3, the bonding cohesion reduces significantly in the 55°–70° contact planes, but does not change very much in the other planes, which agrees with the evolution plotted in Figs. 11(a) and (b). This also agrees with the shear strain distribution in Figs. 13(2-c) and 13(3-c), because the damage of the shear bonding cohesion is controlled by the amount of plastic shear strain.

In the present model, stresses and bonds in each plane are considered as internal state variables and their different evolution in each individual plane introduces in a natural manner the stress-induced anisotropy. Since many soil properties are stress-dependent (e.g., the plastic hardening parameter $k_p$ and the shear modulus $k_r$ are related to $k_n$ which is stress-dependent), the induced anisotropic behavior during undrained shearing is linked to the stress distribution in each plane. The analyses of local stress-strain behaviors also indicate that for real sample, the debonding process is more located in the direction of shear band.
<figure><img src="image_18.png" /><figcaption>**Fig. 13** Local stresses, strains and bonds distributions in rose diagram. (a) For normal stress; (b) for shear stress; (c) for shear strain; (d) for shear bond cohesion</figcaption></figure>

# 6 Conclusions

Cementation among clay clusters is formed when a certain amount of cement is mixed with saturated clay. The influence of the cement content on the mechanical properties of cement-treated clays was investigated. The debonding effects on the stress-strain behavior of cement-treated clay were also discussed. Based on experimental

results, a micromechanical stress-strain model was developed considering explicitly the cementation at inter-cluster contacts. The inter-cluster bonding is introduced by considering two aspects: an additional cohesion in shear sliding and a higher yield stress in normal compression.

Isotropic compression and undrained triaxial tests under various confining stresses on cement-treated Ariake clay and Singapore clay with different cement contents were
simulated to evaluate the capability of the present model to take into account the changes in mechanical properties linked to the cement content. Comparisons between experimental results and simulations demonstrate the model's capability of reproducing correctly the bonding (cementations) and debonding (destructuring) behavior of cement-treated clay with induced anisotropy.

The predicted behavior in contact planes was examined in the case of undrained triaxial tests with different cement contents. The local stress-strain responses in contact planes showed that every contact plane was mobilized to a different degree. A few active contact planes contributed mainly to the deformation of the assembly, while most contact planes contributed only slightly to the specimen straining. Therefore, the local strains are highly non-uniform. It has been shown from the rose diagrams that the shape of the contact stresses, strains and bond distributions, controlling the evolution of local soil properties, changes throughout the mechanical loading, which clearly indicates the development of anisotropy induced by the externally applied loads.

**Acknowledgements** This study was financially supported by the National Natural Science Foundation of China (Grant Nos. 41240024 and 50879058), the Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20110073120012), the Shanghai Pujiang Talent Plan (Grant No. 11PJ1405700), and the ANR-NSFC project RISMOGEO (Grant No. 51161130523).

## Appendix—Micromechanical model of clay

Yin et al. [33] assumed that clay material can be considered as an assembly of clusters. The deformation of a representative volume of clay is generated by mobilizing and compressing of all clusters. Thus, the stress-strain relationship can be derived as an average of the deformation of all local contact planes.

For the $a^{th}$ contact plane, the local forces $f_f^a$ and the local movements $\delta_i^a$ can be denoted as follows: $f_f^a = \{f_n^a, f_s^a, f_t^a\}$ and $\delta_i^a = \{\delta_n^a, \delta_s^a, \delta_t^a\}$ (see Fig. 3 for local coordinate system). To obtain a more direct comparison between the local behavior and the overall stress-strain behavior, we define a local normal stress $\sigma^a = f_n^a Nl/3V$ and a local shear stress $\tau^a = f_r^a Nl/3V$, where $l$ is branch length and $N/V$ is the total number of contact per unit volume. The corresponding local normal strain is defined as $\epsilon^a = \delta_n^a/l$ and a local shear strain is defined as $\gamma^a = \delta_r^a/l$.

### 1) Inter-cluster behavior

For clay, the **elastic stiffness** between two clusters is given by

$$k_n = k_{n0} \left(\frac{\sigma}{\sigma_{at}}\right)^n; k_r = k_{rR} k_{n0} \left(\frac{\sigma}{\sigma_{at}}\right)^n, \quad (A1)$$

where $\sigma_{at}$ is the reference value equal to the atmosphere pressure 101.325 kPa, $\sigma$ is the contact stress in normal

direction. $k_{n0}, k_{rR}$ and $n$ are material constants.

For **friction sliding**, the yield function $F_1$ is assumed to be of Mohr-Coulomb type, defined in a contact-stress space (e.g., $\sigma, \tau_s, \tau_t$), as follows:

$$F_1(\sigma, \tau, H_1) = \frac{\tau}{\sigma} - H_1(\gamma^p), \quad (A2)$$

where $H_1(\gamma^p)$ is a hardening/softening function. Note that the shear stress $\tau$ and the rate of plastic sliding $\gamma^p$ can be defined as

$$\tau = \sqrt{\tau_s^2 + \tau_t^2} \quad \text{and} \quad \gamma^p = \sqrt{(\gamma_s^p)^2 + (\gamma_t^p)^2}. \quad (A3)$$

The hardening function is defined by a hyperbolic curve in $H_1 - \gamma^p$ plane, which involves two material constants: $\phi_p$ and $k_{p0}$.

$$H_1 = \frac{k_{p0} \tan \phi_p \gamma^p}{\sigma \tan \phi_p + k_{p0} \gamma^p}. \quad (A4)$$

The initial slope of the hyperbolic curve, $k_{p0}$, is assumed to relate $k_n$ by a constant $k_{pR}$:

$$k_{p0} = k_{pR} k_n. \quad (A5)$$

The **shear-induced dilation/contraction** is postulated as follows:

$$\frac{d\epsilon^p}{d\gamma^p} = D_a \left|\tan \phi_\mu - \tan \phi_m\right| \left(\frac{\tan \phi_m}{\tan \phi_\mu}\right)^{D_b} \left(\frac{e}{e_c} - 1\right), \quad (A6)$$

where $\tan \phi_m = \tau/\sigma; D_a$ and $D_b$ are material constants for stress dilatancy ($D_b = 1$ can be generally assumed). In Eq. (A6), the $\tan \phi_\mu$ is internal friction angle.

The second plastic law for **normal compression** of inter-cluster contact is expressed as follows:

$$F_2(\sigma, \kappa_2) = \sigma - H_2(\epsilon^p) \quad \text{for} \quad \sigma > \sigma_{pc}, \quad (A7)$$

where $\sigma_{pc}$ is compression yield stress. The normal plastic strain $\epsilon^p$ is generated from the 2nd yield surface only. The hardening function $H_2(\epsilon^p)$ is defined as

$$H_2 = \sigma_{pc0} \exp \left( \frac{\epsilon^p}{c_p} \right), \quad (A8)$$

where $c_p$ is the compression index for the compression curve plotted on $\epsilon^p - \log\sigma$ plane.

The final incremental stress-strain relationship of the inter-cluster contact can be derived, given by

$$\dot{\sigma}_i^a = k_{ij}^{ap} \dot{\epsilon}_j^a \quad \text{or} \quad \sigma = K\epsilon. \quad (A9)$$

### 2) Influence of density state

A density state $\varphi = e_c/e$ is introduced in reference to the void ratio at critical state. The void ratio corresponding to this state is $e_c$. The Critical State Line can be written as follows for clay
$$e_c = e_{cr0} - \lambda \ln \left( \frac{p'}{p_{\text{ref}}} \right), \quad (A10)$$

where $p'$ is the mean effective stress of the packing; $\lambda$ is the slope of critical state line in e-logp' plane, $e_{cr0}$ is the reference void ratio at critical state corresponding to $p_{\text{ref}} = 100$ kPa.

The resistance can be related to the state of packing void ratio e by

$$\tan \phi_p = \left( \frac{e_c}{e} \right)^m \tan \phi_\mu, \quad (A11)$$

where m is a material constant (m = 1 can be generally assumed).

### 3) Overall stress-strain relationship

During the integration process, the static hypotheses are used to obtain the relation between the strain of assembly and inter-cluster strain

$$\dot{u}_{j,i} = \sum_{\alpha=1}^{N} \dot{\gamma}_{j}^{\alpha} n_{k}^{\alpha} B_{ik}^{\alpha}, \quad (A12)$$

where $\dot{\gamma}_j$ is the local strain between two contact clusters; $n_k$ is the unit vector of the branch joining the centers of two contact clusters, and N is the total number of contacts, over which the summation is carried out. The tensor $B_{ik}^\alpha$ in Eq. (A12) is defined as $B_{ik}^\alpha = A_{ik}^{-1} (I^\alpha)^2$, where the fabric tensor

$$A_{ik} = \sum_{\alpha=1}^{N} l_{i}^{\alpha} l_{k}^{\alpha}. \quad (A13)$$

Using the principle of energy balance, which states the work done in a representative volume element equal to the work done on all inter-cluster planes within the element

$$\sigma_{ij} \dot{u}_{j,i} = \frac{1}{V} \sum_{\alpha=1}^{N} f_{j}^{\alpha} \dot{\delta}_{j}^{\alpha} = \frac{3}{N} \sum_{\alpha=1}^{N} \tau_{j}^{\alpha} \dot{\gamma}_{j}^{\alpha}; \quad (A14)$$

and using Eq. (A12), the local stress on the $\alpha$th contact plane is derived as follows:

$$\dot{\tau}_{j}^{\alpha} = \frac{N}{3} \dot{\sigma}_{ij} B_{ik}^{\alpha} n_{k}^{\alpha}. \quad (A15)$$

The stress increment $\dot{\sigma}_{ij}$ can be obtained by the contact forces and branch vectors for all contacts. In terms of local stress, it is

$$\dot{\sigma}_{ij} = \frac{1}{V} \sum_{\alpha=1}^{N} f_{j}^{\alpha} l_{i}^{\alpha} = \frac{3}{N} \sum_{\alpha=1}^{N} \tau_{j}^{\alpha} n_{i}^{\alpha}. \quad (A16)$$

Using Eqs. (A9), (A12), and (A15), the following relationship between stress and strain can be obtained:

$$\dot{u}_{i,j} = C_{ijmp} \dot{\sigma}_{mp}, \quad (A17)$$

$$\text{where } C_{ijmp} = \frac{N}{3} \sum_{\alpha=1}^{N} (\bar{k}_{jp}^{ep})^{-1} n_{k}^{\alpha} n_{n}^{\alpha} B_{ik}^{\alpha} B_{mn}^{\alpha}. \quad (A18)$$

A numerical calculation with an iterative process is used to carry out the summation in Eq. (A18).

## References

1. Tatsuoka F, Kobayashi A. Triaxial strength characteristics of cement-treated clay. In: Proceedings of the 8th European Conference on SMFE. Helsinki, 1983, 8(1), 421–426

2. Tremblay H, Leroueil S, Locat J. Mechanical improvement and vertical yield stress prediction of clayey soils from eastern Canada treated with lime or cement. Canadian Geotechnical Journal, 2001, 38(3): 567–579

3. Horpibulsuk S, Miura N, Bergado D T. Undrained shear behavior of cement admixed clay at high water content. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(10): 1096–1105

4. Kamruzzaman A H M, Chew S H, Lee F H. Structuration and destructuration behavior of cement-treated Singapore marine clay. Journal of Geotechnical and Geoenvironmental Engineering, 2009, 135(4): 573–589

5. Mitchell J K. Soil improvement—State of the art report. In: Proceedings of the 10th International Conference on Soil Mechanics and Foundation Engineering. Balkema, The Netherlands, 1981, 4: 509–565

6. Çokça E. Use of Class C flyashes for the stabilization of an expansive soil. Journal of Geotechnical and Geoenvironmental Engineering, 2001, 127(7): 568–573

7. Yin Z Y, Hattab M, Hicher P Y. Multiscale modeling of a sensitive marine clay. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(15): 1682–1702

8. Yin Z Y, Chang C S, Hicher P Y, Wang J H. Micromechanical analysis for the behavior of stiff clay. Acta Mechanica Sinica, 2011, 27(6): 1013–1022

9. Mitchell J K. Fundamentals of Soil Behaviour. 2nd ed. New York: Wiley, 1992

10. Tatsuoka F, Uchida K, Imai K, Ouchi T, Kohata Y. Properties of cement treated soil in Trans-Tokyo Bay highway project. Ground Improvement, 1997, 1(1): 37–57

11. Balasubramaniam A S, Lin D G, Sharma S S, Kamruzzaman A H M, Uddin K, Bergado D T. Behavior of soft Bangkok clay treated with additives. In: Proceedings of the 11th Asian Regional Conference on Soil Mechanics and Geotechnical Engineering. Balkema, The Netherlands, 1999, 11–14

12. Burland J B. On the compressibility and shear strength of natural soils. Geotechnique, 1990, 40(3): 329–378

13. Leroueil S, Vaughan P R. The general and congruent effects of structure in natural soils and weak rock. Geotechnique, 1990, 40(3): 467–488

14. Liu M D, Carter J P. Virgin compression of structured soils. Geotechnique, 1999, 49(1): 43–57

15. Lee K, Chan D, Lam K. Constitutive model for cement treated clay in a critical state framework. Soil and Foundation, 2004, 44(3): 69–77

16. Horpibulsuk S, Liu M D, Liyanapathirana D S, Suebsuk J. Behaviour of cemented clay simulated via the theoretical framework
of the structured cam clay model. Computers and Geotechnics,
2010, 37(1-2): 1-9

17. Batdorf S B, Budianski, B A. Mathematical theory of plasticity based on concept of slip. NACA Technical Note TN, 1949, 1871

18. Calladine C R. Microstructural view of the mechanical properties of natural clays. Geotechnique, 1971, 21(4): 391-415

19. Rothenburg L, Selvadurai A P S. Micromechanical definitions of the Cauchy stress tensor for particular media. In: Selvadurai APS, ed. Mechanics of Structured Media. Amsterdam: Elsevier, 1981, 469-486

20. Pande G N, Sharma K G. Multi-laminate model of clays—a numerical evaluation of the influence of rotation of the principal stress axis. In: Desai C S, Saxena S K, eds. Symposium on Implementation of Computer Procedures and Stress-Strain Laws in Geotechnical Engineering. Chicago, Durham, NC: Acorn Press, 1982, 575-590

21. Jenkins J T. Volume change in small strain axisymmetric deformations of a granular material. In: Satake M, Jenkins J T, eds. Micromechanics of Granular Materials, Amsterdam: Elsevier, 1988, 143-152

22. Chang C S, Liao C. Constitutive relations for particulate medium with the effect of particle rotation. International Journal of Solids and Structures, 1990, 26(4): 437-453

23. Bazant Z P, Xiang Y, Ozbolt J. Nonlocal microplane model for damage due to cracking. Proceedings of Engineering Mechanics, 1995, 2: 694-697

24. Chang C S, Gao J. Second-gradient constitutive theory for granular material with random packing structure. International Journal of Solids and Structures, 1995, 32(16): 2279-2293

25. Emeriault F, Cambou B. Micromechanical modelling of anisotropic non-linear elasticity of granular medium. International Journal of Solids and Structures, 1996, 33(18): 2591-2607

26. Nicot F, Darve F. A multiscale approach to granular materials. Mechanics of Materials, 2005, 37(9): 980-1006

27. Chang C S, Hicher P Y. An elastoplastic model for granular materials with microstructural consideration. International Journal of Solids and Structures, 2005, 42(14): 4258-4277

28. Cudny M, Vermeer P A. On the modelling of anisotropy and destructuration of soft clays within the multi-laminate framework. Computers and Geotechnics, 2004, 31(1): 1-22

29. Galavi V, Schweiger H. A multilaminate model with destructuration considering anisotropic strength and anisotropic bonding. Soil and Foundation, 2009, 49(3): 341-353

30. Chang C S, Hicher P Y, Yin Z Y, Kong L R. An elasto-plastic model for clay with microstructural consideration. Journal of Engineering Mechanics, 2009, 135(9): 917-931

31. Yin Z Y, Chang C S. Microstructural modelling of stress-dependent behaviour of clay. International Journal of Solids and Structures, 2009, 46(6): 1373-1388

32. Yin Z Y, Chang C S. Non-uniqueness of critical state line in compression and extension conditions. International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(10): 1315-1338

33. Yin Z Y, Chang C S, Hicher P Y, Karstunen M. Micromechanical analysis of kinematic hardening in natural clay. International Journal of Plasticity, 2009, 25(8): 1413-1435

34. Lagioia R, Nova R. An experimental and theoretical study of the behaviour of a calcarenite in triaxial compression. Geotechnique, 1995, 45(4): 633-648