Numerical and experimental modeling of geomechanical behavior of partially saturated soils

Abstract

The Barcelona Basic Model (BBM) has been implemented in a finite difference-based computer program to simulate the behavior of unsaturated soils subjected to wetting. The BBM implementation was verified using analytical solutions, and the proposed model has been used to simulate the response of a compacted embankment under complete saturation and suction induced conditions. Numerical analyses indicate that considerable amount of total and differential settlements could develop at the top surface of the embankment. BBM is executed into FLAC^2D extending a defined module for modified Cam Clay (MCC) and has been set up an analytical solution for suction-dependent stress and strain. Evaluating the effect of anisotropy and nonlinear apparent tensile strength in unsaturated soils, a modification to BBM formulation has been proposed and optimized by developing numerical analyses to reduce the size of elastic region of loading collapse (LC) curve. Then, an experimental study in the literature is investigated by utilizing comparative curves from BBM and modified BBM indicating well agreement with natural circumstances. As a result of the work presented in this research, finite difference codes with BBM and modified BBM has the capability of simulating the real behavior and is operational being applied to problems associated with earthen structures in unsaturated or partially saturated of expansive soils as a three-phases medium.

Introduction

Constitutive modeling of unsaturated soils has been a main subject of many recent studies of researchers [2, 21, 24, 28, 29, 31, 33, 35, 36, 38]. However, those sophisticated constitutive models are seldom implemented in computer programs for practical applications because of model complexity. The mechanical behavior of soil is dependent on the both current state of stress type of soil, stress history and stress paths. The Barcelona Basic Model (BBM) is a geomechanical constitutive model for obtaining the elasto-plastic behavior of unsaturated soils [8, 9, 17]. The model was first developed and presented in the early1990s as a spread of the modified Cam-Clay (MCC) model to unsaturated soil conditions [2].

In fact, first integrated constitutive model capable to predict the various aspects of partially saturated soil behavior, called loading-collapse (LC) model, was presented by Alonso et al. [2] in which major focus was on volume change behavior, including collapse. The model can describe several typical features of unsaturated-soil mechanical behavior, including wetting-induced swelling or collapse strains, depending on the magnitude of applied stress, as well as the increase in shear strength and apparent pre-consolidation stress with suction [5, 15]. To authenticate the BBM, Wheeler and Sivakumar [34] conducted a series of experimental research on compacted kaolin soil indicating a good agreement with BBM predictions.

In this research, implementation of Barcelona Basic Model into the FLAC is presented analyzing geo-mechanical behavior of unsaturated soil. In the present approach, an existing module for modified Cam-Clay (MCC) model in FLAC^2D is extended and a function as a computational solution for suction-dependent strain and net stress (i.e. total stress minus air pressure: $\sigma_{t} - \sigma_{a}$) in unsaturated soil medium. Results of centrifuge model [39] was used to simulate the suction response of a compacted model as an embankment, evaluating displacements in the model. A three stages procedure is mainly used to provide unsaturated soil properties of different applied matric suctions and different paths by using single soil specimen. Stress paths are produced for loading collapse curves, using two cases of considering predicted BBM-Alonso and experimental results. Then, the effect of anisotropy and none linear apparent tensile strength are investigated. In the next step, a modification to a part of the constitutive model equation will be proposed to control elastic zone in the model in higher suctions. As a result of implementing this modified model in FLAC, two stress paths for LC curve are plotted, comparing with the results from both Alonso model implemented in FLAC and experimental researches in this study. In continuation of this research, an experimental study performed by Cui and Delage [4] in a controlled suction triaxial apparatus is investigated. The stress–strain and volume change plotted from the experimental study, BBM and modified BBM are evaluated and compared. [27] modified the BBM model with respect to over consolidated ratio.

Barcelona Basic Model Formulations

Barcelona Basic Model (BBM) is formulated in triaxial stress space using modified Cam Clay model (MCCM) as the underlying model for saturated condition. Under isotropic stress states the relevant stress variables in the BBM are the mean net stress $\overline{p}$ and matrix suction s defined as follows [23]:

$$\overline{p} = p - u_{a} \quad \quad s = u_{a} - u_{w}$$

(1)

where $\overline{p}$ is the mean total stress, u_a is the pore air pressure and ${\text{u}}_{{\text{w}}}$ is the pore water pressure. In the model, elastic changes of specific volume v are related to changes of $\overline{p}$ and s as

$${\text{d}}\nu_{e} = - \kappa \frac{{{\text{d}}\overline{P}}}{{\overline{P}}} - \kappa_{s} \frac{{{\text{d}}s}}{{s + P_{at} }}$$

(2)

where ${\upkappa }$ and ${\upkappa }_{{\text{s}}}$ are two elastic constants and p_at is atmospheric pressure. For isotropic stress states the onset of plastic strains coincides with a loading collapse (LC) yield curve in the $\overline{p} - s$ plane. Post-yield compression curves at different values of suction define a series of isotropic normal compression lines in the v-$p$ plane. It can be stated that the relationship for normal compression line (NCL) at any constant suction s, is given by the following equation:

$$v = N\left( s \right) - \lambda \left( s \right)ln \overline{p}$$

(3)

where the intercept $N\left( s \right)$ and gradient $\lambda \left( s \right)$ are both functions of suction (see Fig. 1a). Note that the intercept $N\left( s \right)$ in Eq. 3 is defined at $\overline{p} = 1$ kPa and is not the same as the intercept.

Fig. 1

a Post-yield compression curves at different values of suctions, b the shape of the LC yield curve in the $\overline{p }-s \mathrm{plane}$ [33]

defined by Alonso et al. [2] and Josa [21] at a reference pressure p_c. By combining Eqs. 2 and 3, the shape of the LC yield curve in plane, and the development of this shape during yield curve expansion, and pre-consolidation pressure p_c are derived as

$$\frac{{p_{0} }}{{p_{c} }} = \left( {\frac{{p_{0}^{*} }}{{p_{c} }}} \right)^{{\left( {\frac{\lambda \left( 0 \right) - \kappa }{{\lambda \left( s \right) - \kappa }}} \right)}} \quad \quad \frac{{\overline{p}_{0} \left( {\left( s \right)} \right)}}{{p_{c} }} = \left( {\frac{{\overline{p}_{0} \left( 0 \right)}}{{p_{c} }}} \right)^{{\left( {\frac{\lambda \left( 0 \right) - \kappa }{{\lambda \left( s \right) - k}}} \right)}}$$

(4)

where $\overline{p}_{0}$(s) is the isotropic yield value of $\overline{p}$ at a general suction s (see Fig. 1b), $\overline{p}_{0}$(0) is the value of $\overline{p}_{0}$(s) at zero suction, the intercept of the yield curve with the $\overline{p}$ axis, defining the current size of the curve) and p_c is an additional soil constant equals to $\overline{p}_{0}$(0). A major assumption of Alonso et al. [2] and Gens and Potts [14] in deriving Eq. 5 was that the developing shape of the yield curve as it expanded could be traced back to a situation where the yield curve would be a straight vertical line in the s:p plane at the reference stress p_c (see Fig. 1b). $p_{0}^{*}$ is the pre-consolidation pressure at zero suction and as $\overline{p}_{0}$(0) is one of the hardening parameters. Figure 1b illustrates that the LC yield curve is vertical in the s:p plane when $\overline{p}_{0}$(0) = p_c, the yield curve then becomes increasingly inclined as it expands (i.e. as $\overline{p}_{0}$(0) increases). Illogical yield curve shapes are suggested for values of $\overline{p}_{0}$(0) less than the reference or pre-consolidation pressure p_c, and it is therefore important that a low value is selected for the soil constant p_c (significantly lower than the lowest value of $\overline{p}_{0}$(0) likely to occur in a given application). The locations of the normal compression lines for different values of suction (based on Eq. 3) are linked to the developing form of the LC yield curve as it expands based Eq. 5, and it can be shown that the intercept $N\left( s \right)$ of the normal compression line at a given value of suction is:

$$N\left( s \right) = N\left( 0 \right) - \left( {\lambda \left( 0 \right) - \lambda \left( s \right)} \right){\text{ln}}p^{c} - \kappa_{s} {\text{ln}}\left( {\frac{{s + p_{atm} }}{{p_{atm} }}} \right)$$

(5)

where $N\left( 0 \right)$ is the corresponding intercept at zero suction. The seven soils constants required to model behavior under isotropic stress states with the BBM are therefore κ, κ_s, $\lambda \left( 0 \right)$, $N\left( 0 \right)$, r, β and p_c. Three further constants are required to complete the model for triaxial stress states.

The extension of the modified Cam Clay (MCC) model to unsaturated media and non-isothermal conditions adds to the complexity and the number of required input parameters. Five parameters are needed to define the MCC model, and nine additional parameters are necessary for the BBM model. In the FLAC implementation of the BBM, it was also developed some FISH routines, such as one for calculating suction strain. A function for suction alternations were programmed with FISH codes which simulate capability of yield functions in three-phase unsaturated regions

$$\lambda \left( s \right) = \lambda \left( 0 \right)\left[ {\left( {1 - r} \right)\exp \left( { - \beta .s} \right) + r} \right]$$

(6)

This equation explains suction increase of yield function. $\lambda \left( 0 \right)$, r and β are soil constants and key parameters controlling virgin loading under isotropic stress state. Alonso et al. [2] Gens and Potts [14] assumed that λ(s) decreases with increasing suction as displays in Fig. 2, therefore the value of r was less than 1.

Fig. 2

Compression curves: a In the p-υ plane; b In the s-υ plane [2]

The BBM includes two independent stress variables: the net stress (σ − u_a) and matric suction (u − u_w), in which σ is the total stress, and u_a and u_w are the pore air and pore-water pressures, respectively. The following assumptions were made for the implementation of BBM: (1) pore-air pressure is zero relative to atmospheric pore-air pressure, which is true for most natural ground conditions [38], therefore, net stress is equivalent to the total stress,and (2) matric suction is a parameter that influences the stiffness and strength of unsaturated soils during the course of wetting. The implementation of BBM in FLAC subsequently is described briefly, and more details are reported by Zheng et al. [39]. In the BBM, both the net stress and matric suction can induce strain in the soil [25]. Strains induced by net stress (in the p-υ plane) and suction (in the s-v plane) are considered separately, as shown in Fig. 2. The incremental volumetric strains induced by net stress are divided into two components: elastic strain $d\varepsilon_{v}^{e}$, and plastic strain $d\varepsilon_{v}^{p}$,as reported by Alonso et al. [2]

For wetting-induced deformation problems, the wetting path is a suction-unloading process [19]. At any point during wetting, the current suction s is equal to or smaller than the maximum suction s₀ ever experienced by the soil in the past. Therefore, the soil remains in the elastic region in the s-v plane during the course of wetting. Recent researches compute the suction-induced strain by adding an equivalent mean stress increment to the net stresses, which is also adopted in this study and suits for the implementation in FLAC.

Verification

Numerical simulations were conducted for self-weight, complete saturation and suction variations condition, and the implemented of BBM in FLAC has been verified using analytical solutions in the literature. The numerical results comprising model surface and centerlines’ settlements are presented. BBM parameters for the numerical simulation were selected from Alonso’s research which are summarized in Table 1. It should be mentioned that the suction-induced routes due to wetting process are simulated by suction decrement from a determined value to the minimum amount (i.e. 200 kPa to 0). This numerical analysis can simulate the reduction changes in suction in a relatively simple way. Figure 3 demonstrates stress paths by volumetric deformation due to wetting process from Alonso’s experimental study, case 1. Comparative curves for different stress paths indicated in Fig. 4. The change of specific volume v by various net mean stresses has been plotted. The results from Alonso et al. [2] and numerical analysis in the present study specified regular lines and dash lines respectively. The path starts with ABFH line with swelling inside the elastic zone. The sample experienced compaction in paths AEFH and AEGH. Comparison of the results shown in Fig. 4, indicate a good agreement between the predicted [2] and analytical results with FLAC with respect to both magnitudes and trends.

Table 1 BBM parameters as adapted from Alonso et al. [2]

Fig. 3

Stress path in literature case study [2]

Fig. 4

Comparison of FLAC prediction and BBM analytical results in the stress paths

Numerical model of an unsaturated embankment

In order to evaluate the proposed model, an embankment in different stages of construction with fully saturation is considered [11]. BBM parameters for the reference soil reported by Alonso et al. [2] are summarized in Table 1.

The stress paths discussed earlier are most commonly used to determine parameters which are adapted from some isotopically triaxial tests proposed by Alonso et al. [2]. In this section, a compacted silty clay from [2, 21] with BBM parameters was selected as the model properties. This type of soil is categorized as symbol group of ML (silt) in ASTM D-2487. Figure 5 shows the geometry and finite difference grid for embankment model.

Fig. 5

Geometry and finite difference grid for the embankment model

The behavior of foundation saturated soil has been considered as Mohr–Coulumb model and was assumed to be in a saturated condition. The height of the model was 10 m with 1:2 slope and depth of soil as a foundation was 5 m.

The embankment can be partially modeled because of its symmetry. The model centerline was placed at the embankment crown along the boundary-line. The embankment slopes are fixed for lateral displacement. However, in the bottom side both displacements in x and y directions are fixed. The model was simulated utilizing the BBM relations with an initial suction amount of 100 kPa. Garcia [13] and Hatami et al. [16] investigated on compacted silty soil and conclude a comparative suction amount. In the first step of the model construction the foundation are reached equilibrium under its weight because of gravitational stress. As real circumstances, the embankment was constructed in 10 steps, considering initial suction of 100 kPa. And while, suction value alters from amount of 100 kPa to zero, the whole embankment is subjected to impermanent equilibrium. In other words, the 10 m embankment is simulated in 10 stages to depict the stages of construction correctly.

Numerical analyses

For verifying BBM-Alonso formulations in FLAC, analytical solutions have been utilized with those in the same conditions in literature. For the better investigation of BBM in FLAC, settlements, displacements in different cases are evaluated to verify this study in comparison with literature. Figure 6 depicts a comparative scatter for settlements of embankment surface at different stages. In Zheng et al. [39] study, settlement of the model surface is relatively 60 mm, and in the present study, this value is 53 mm. After complete saturation, the settlement of embankment surface is approximately two or three times larger than that the last situation, amount of differential settlement, reached approximately 239 mm and 250 mm, in this study and the literature, respectively, as shown in Fig. 6 by the scatter in different distance from the centerline. Hence, the amounts of suction-induced settlement is nearly 190 mm. Figure 6 represents the settlement along the longitudinal stages of the embankment in two conditions, self-weight loading and full saturation as mentioned above. The scatters illustrate that y-displacement of embankment at different cases from the present study are in a close agreement with a selected published research [1]

Fig. 6

Comparison of settlements of embankment surface at different distances from centerline

It can also be seen in Fig. 7 showing settlements of embankment centerline along elevation at different stages. The maximum settlement regarding to self-weight occurred at mid-height of the embankment as the model was constructed gradually and the accumulative displacement was reported since the placement of each level. The results in both studies [37] with the present research) match very well, as are being observed in Fig. 7.

Fig. 7

Comparison of settlements in different stages in three cases of self-weight loading, full saturation and wetting-induced

A modification to Barcelona Basic Model relations

Stress paths

The first comprehensive constitutive model which can be defined for predicting different aspects of unsaturated or partially saturated soil, namely loading-collapse (LC) model in which its main focus was on specific volume behavior, containing collapse [10]. A new model proposed by Alonso et al. [2] resulted from quantitative version of the LC model, is Barcelona Basic Model (BBM), where changes in specific volume and strength behavior were anticipated within the framework of elasto-plastic theory. This model was based on the experimental model over compacted kaolin. Alonso et al. [2] calibrated their tests according to Kubo and Karube [22] and Josa et al. [20] experimental studies. The soil properties for their tests were: PI = 12% PL = 27% LL = 39% with largest grain size of 0.04 mm. The stress paths and initial position of the both yield curves loading collapse and suction increase in BBM codes are depicted in Fig. 8a provided by Alonso et al. [2], case 1. Path ABFH occured inside the elastic region and persuaded a little swelling. In stress path AEFH the sample was compacted to a specific value for net mean stress, as indicated in Fig. 8a. In this path, specific volume decreased from 1.82 to 1.76 and similarly path AEGH started at a confining pressure with amount of 600 kPa and yielded a specific volume reduction from 1.79 to 1.73. The sample soil met low value of swelling for smaller stresses, in the elastic medium. On the other hand, after yielding at bigger stress planes, the curve fell down.

Fig. 8

a Response of model to suction-Specific volume-net mean stress. b. Response of model to suction- Stress paths

In all three cases, loading–collapse is induced by the movements of the LC curve from its initial position LCi to its final position LC_H. FLAC results in Fig. 8b agree with analytical results presented in Fig. 11 from Alonso et al. [2].

This procedure is from the result of triaxial test with matric suction control. It involves imposing and controlling positive air and water pressure in the soil sample (higher than atmospheric pressure) and allows controlling the matric suction in the soil [6]. It is worthy to mention that the behavior of unsaturated soils is related by the soil suction, ψ. Soil suction is composed of matric suction and osmotic suction. Matric suction is defined as the difference between pore air pressure, u_a, and pore water pressure, u_w.

$$\psi = u_{a - } u_{w}$$

(7)

Determination of matric suction distribution in unsaturated medium represents a major step, necessary for further geotechnical applications regarding specific unsaturated soil behavior. A suction control triaxial apparatus can be used to perform test under unsaturated condition by a special technique for controlling pore water and air pressure in order to preserve matric suction within the specimen. Some control devices and measurement instruments are equipped in the apparatus, such as: external and internal load cell to measure deviatoric stress; pressure transducer to measure cell pressure, pore water pressure transducer and volume change control [32]. In fact this set is a modified version of an ordinary triaxial apparatus for testing on unsaturated soils that is shown in Fig. 9.

Fig. 9

Schematic of a modified triaxial apparatus for testing unsaturated soils [6]

As mentioned before, Alonso provided his model based on a fine-grained soil with a maximum nominal aggregate size of 0.04 mm. The soil type used for the test is a sedimentary which is abundant in the eastern parts of North America. The amount of natural moisture in the soil varies between 20 and 29%. In-situ porosity and specific volume of this low plasticity soil are respectively 1.39 and 16 kN/m².

The variables in the BBM are the mean net stress $\overline{p}$ and matric suction s, are defined as $\overline{p}$ is the mean total stress, u_a is the pore air pressure and $u_{w}$ is the pore water pressure.

A three stages procedure is mainly used to provide unsaturated soil properties of different applied matric suction and various paths by using single soil specimen. BBM parameters suggested for this soil is given in Table 2 and a typical set of series for each stage is presented in Table 3.

Table 2 BBM parameters selected for the soil specimen [2]

Table 3 Schematic stresses and matric suction for different stages in triaxial stress test

According to the available laboratory results, the stress path ABFH and AEFH were investigated (Fig. 10). The difference between the above stress paths is only in the mean stress where the soil is saturated, as mentioned before. An issue to be considered is porosity and initial specific volume of the soil

$$e = 1.44 \Rightarrow \nu = 1 + e = 2.44$$

(8)

with the new soil parameters summarized in Table 1. Implementation of BBM formulation in FLAC, and the new loading collapse (LC) curve are plotted in Fig. 11.

Fig. 10

Stress paths for the present soil type

Fig. 11

Comparison between measured and predicted model results for stress path AA₁BH

Figures 11 and 12 illustrate LC diagram from the model with BBM parameters of silty sand with low plasticity (LCT1) showing that the change in the volumes in the model before the saturation phase of the soil is smaller than the change in the volumes obtained from the experiment. However, after saturation, the model and the experimental results are in good agreement. The reason is that BB-model converts to modified Cam-Clay when soil is fully saturated. This issue is also observed in the comparative graph operated by Alonso and Josa in 1990 (Fig. 13). One reason could be that the deformation coefficient is larger for the mean stress. Another reason may be due to the large elastic area proposed by the model. However, there seems to be a need for modification of the model to obtain more realistic results which match with in-situ soils. In fact there are high suction and anisotropy in natural soils.

Fig. 12

Comparison between measured and predicted model results for stress path AEFH

Fig. 13

Comparison between measured and predicted results for partially unsaturated clay [2]

Nonlinear apparent tensile strength approach

A distinguished feature of the original BBM is the assumption of a linear apparent tensile strength (ATS) approach in p: s plane, as a full set of apparent tensile strength values P_s(s) for each associated value of suction; where, k = linear rate of increase of P_s(s) with suction. The exact location of yield surface in a partially saturated soil subjected to suction controlled shearing is depended on the final position of critical state lines, shear strain increment and variation apparent tensile strength by suction [6]. A plenty of constitutive models have been developed to predict the shear strength of these types of soils. These constitutive models including BBM described a linear increase of apparent tensile strength (ATS) with suction. However, this has been found unreasonable and suitable only for low values of suction. At high value of suction, while the soil become completely dry the increase of apparent tensile strength due to suction leads to zero [7]. This issue indicates that the apparent tensile strength commonly grows with suction and once reaches a maximum value, it decreases regularly tending to a low value at very high value of suction. Some experimental researches present the effect of a non-linear apparent tensile strength with suction (ATS) for a wide range of soils. Hoyos et al. [18] performed a series of formal triaxial compression tests on clayey sand specimens under different constant-suction states. Experimental loading- collapse approach determined on p₀–s plane was well in confirmed with the one obtained using LC equation. However, the acquired ATSL was nonlinear. Hoyos et al. [18] expressed relationship for function p_s with suction s as, $p_{s} = - ks^{m}$ in which, m and k control the rate of change of $p_{s}$ by curvature and ATSL suction, respectively. The best proper amount in experimental data was applied utilizing the method of least squares to assess the model parameters k and m in the above equation.

Anisotropy

A set of triaxial shear tests were performed on compacted silt commencing with consolidation phase followed by shearing with various stress ratios ƞ [4]. This experimental research indicates that several soils express anisotropy of mechanical behavior arising from the soil structure anisotropy which may be emanated throughout the formation of the soil. This anisotropy can be changed afterward by plastic straining which affects re-arrangement of the initial compound. An anisotropy elasto-plastic critical state constitutive model S-CLAY 1 for saturated soils was developed considering plastic behavior of anisotropy. This behavior was brought up via an inclined yield surface and a rotational component of hardening to the shape of the development or elimination of the initial anisotropy throughout plastic deformation [33]. An anisotropic elastoplastic model for partially saturated soil through combination of BBM and S-CLAY 1 model as Anisotropic Barcelona Basic Model [30]. Only anisotropy of plastic strains was modelled by Anisotropy Barcelona Basic Model (ABBM) whereas the elastic behavior is isotropic and defined by formulation identical to BBM. Gallipoli et al. [12] proposed some modifications to Anisotropic Basic Model (ABBM) developed by Stropeit et al. [30]. Yield surface for any constant suction is defined in p₀—q plane as:

$$f_{1} = \left( {q - \alpha p^{\prime}} \right)^{2} - \left( {M^{2} - \alpha^{2} } \right)\left( {p^{\prime} + \frac{M}{M - \alpha }p_{s} } \right)\left( {p_{mc} - p^{\prime}} \right) = 0$$

(9)

where $\alpha$ is yield surface inclination which controls the rotational component of hardening as $p_{mc}$ size of yield surface which controls the volumetric component of hardening. Terms M, p_s, q and p₀ are critical state line slope, function of suction, deviatoric stress and initial value of hardening parameter as mentioned in BBM formulation.

A proposed modification for model-elastic zone

From all topics and discussions at this stage It would be concluded that a feasibility study for model modification can be accomplished. The experimental results in the literature show that it seems the elastic region of the model is wide for high amounts of matric suction and it is possible to obtain more realistic results if this area to be reduced. LC yield function were investigated in stress paths, such that loading collapse (LC) yield function can be modified. According to similar LC yield curve expression (Eq. 4) as for BB-Model given by

$$\frac{{p_{mo} }}{{p_{c} }} = \left( {\frac{{p_{m0}^{*} }}{{p_{c} }}} \right)^{{\left( {\frac{\lambda \left( 0 \right) - \kappa }{{\lambda \left( s \right) - \kappa }}} \right)}}$$

(10)

The hardening component $p_{mo}$ which depends on suction, govern the size of constant suction yield surface, where $p_{m0}^{*}$ defines the size of the yield curve at s = 0 as shown in Fig. 14.

Fig. 14

Constant suction cross section of yield surface in BBM

Therefore, a factor should be considered to control the elastic area in high suction and to make a fewer changes in lower ones. For this reason, this coefficient could be inversely related to suction and multiplied as a function of suction in LC equation. If S_max and S₀ are defined as maximum amount for matric suction and initial suction, respectively, m factor can be proposed as

$$m = \frac{{S^{\prime} - S_{0} }}{{S^{\prime} + S_{0} }}$$

(11)

Practically, the value of S_max shouldn’t exceed a certain value because mathematically its value may reach very close to one and loses its performance. To optimize the value of S_max, the implementation of a model with BBM parameters in FLAC is applied. A new simulation analysis of the model was conducted to verify implementation of the modified BBM with m factor. A compacted silty soil with calibrated BBM parameter [2] was selected as fill for the embankment (Table 1) Fig. 14 demonstrates the settlement along the longitudinal stages of the embankment in three cases of BBM Alonso, normally LC equation and modified curve with m factor in the present study. As displayed in Fig. 15, the curve based on the modified relation are more optimized and better tendency to the real actual situation. The most acceptable results were obtained when S_max was approximately 4.5 to 5.5 times the Sʹ which is the highest suction soil experienced as

$$S^{\prime} = \left( {4.5\sim 5.5} \right)S_{max}$$

(12)

Fig. 15

Settlements of embankment surface at different stages, both in this study (ordinary BBM and modified one) and the literature [39]

On the other hand, there are restrictions to the model modifications. This means, in low amount of matric suctions, the simulation results are satisfied, therefore no changes are acceptable in this case. Afterwards, with providing several data from literature and utilizing in numerical computations for none-linear regression [26], m factor approximates to the following formulation

$$m = \left( {\frac{{S^{\prime} - S_{0} }}{{S^{\prime} + S_{0} }}} \right)^{{\left[ {\log \beta r} \right]^{0.5} }}$$

(13)

where β and r are parameters defining the maximum increase of soil stiffness with suction representing the maximum soil stiffness.

As the results of this session, LC equation can change to the following relation:

$$\frac{po}{{p_{c} }} = m\left( {\frac{{p_{0}^{*} }}{{p_{c} }}} \right)^{{\left( {\frac{\lambda \left( 0 \right) - \kappa }{{\lambda \left( s \right) - \kappa }}} \right)}}$$

(14)

Figure 16 illustrates the results of a parametric study performed by Zheng et al. [39] to investigate the influences of suction. Matric suction generally decreases as getting near the groundwater. The initial matric suctions are 100 and 200 kPa. Figure 13a indicates settlements of the embankment surface at different stages. Results indicate that higher amount of settlement would be anticipated for a compacted kaolin placed at wetter conditions (i.e. lower initial suction) under self-weight compression. The results from BBM and MBBM (modified BBM) implemented in FLAC, figured to compare with the previous researches. In both, the embankment with lower initial suction is expected to undergo a greater amount of self-weight settlement. However, this amount is not as greater as these two cases and have been optimized. For S_i = 200 kPa, the surface settlement after full saturation was more than 70 mm larger than that for S_i = 100 kPa. Meanwhile, the differential settlement for MBBM is less than 50 mm. Lower amount of settlements is another aspect from MBBM considering m ratio in BBM formulation. Figure 15b depicts that wetting-induced settlement of the model surface increases moderately, as initial soil suction increases. Considering MBBM as the constitutive model cause the embankment surface to be settled with lower amounts. Results in this figure confirm that for BBM implementation case, the wettest unsaturated embankment with S_i = 100 kPa exhibited the largest amount of settlement whereas the model indicates the smallest amount of settlement in which MBBM considers for constitutive model. Therefore, m ratio is an effective factor for BBM modification in wetting process of unsaturated silty clays for settlement control problems.

Fig. 16

Comparison of initial suction effect on settlements of the model surface in various situations for: a Self-weight settlements; b Wetting-induced settlements

Capabilities of the proposed modified BBM

In this session, the results from experimental study and original model implemented in FLAC, compare with m ratio to demonstrate the performance of modified BBM. The soil is the same with parameters summarized in Table 1. Figures 17 and 18 display comparative LC diagrams in different situations of BBM-Alonso prediction, experimental and modified BBM results.

Fig. 17

Comparison among measured, predicted and modified model results for stress path AA₁BH

Fig. 18

Comparison among measured, predicted and modified model results for stress path ABFH

As specified in the figures, describing m ratio for a LC equation from BBM codes as a modification, has an important effect on model behavior, such that results from modified model are closer to experimental results compared with the original model. The figures depicts that the modified model has proved better performance in predicting volume changes.

Assessment of plastic behavior of an unsaturated compacted silt

To investigate the modified model, an experimental research performed by Cui and Delage [4] in an osmotically controlled suction triaxial apparatus. The experimental behavior of compacted silt by conducting tests on aeolian silt from eastern region of Paris, are studied. In their experimental studies, the volume change by shear Isotropic loading test was focused. Also, pre-consolidation behavior by increasing suction was studied, and a nearly matric suction hardening phenomenon was revealed, together with a non-associated flow rule. In the present study the situation was simulated in finite difference programming in which BBM codes are implemented. The unsaturated model parameters for Jossigny silt utilized in programming are summarized in Table 4. In the laboratory [4] different cell pressures (50, 100, 200, 400 and 600 kPa) and four suction values (200, 400, 800 and 1500 kPa) were generally applied. The comparative stress–strain curves at s₃ = 50 kPa under higher suctions (800 and 1500 kPa) among the experimental study, BBM and modified BBM in FLAC are presented in Fig. 19a. For these amounts of suctions, an increasingly pronounced maximum values is observed at a reducing level of strain, displaying an increasing fragility of the sample with suction. The same experimental study showed that peaks corresponding to the development of shear planes, proportional to failure by strain. At strains larger than 8% a residual condition seems to be reached, with a constant value of deviator stress wisely independent of the suction. For BBM curve, the strain increase process is very slow and fragility occurs before the axial strain of 6%. For higher suction (1500 kPa), the greater deviator stresses, the greater the distance between the actual values and the model. On the other hand, there are well agreements between MBBM and the experimental results. In addition to that, failures took place approximately at the same points as occurred in the experiments. In experimental studies, the volume change curves (Fig. 19b) are affected by axial strains. For higher suctions (800 kPa and 1500 kPa), contraction followed by dilation is observed: the higher the suction, the larger is the dilatancy. Of course, volume change curves are valid only before any peak in the deviator stress curve. All curves are well ordered and grow regularly with increasing suction.

Table 4 BBM parameters selected as adapted from Cui and Delage [4]

Fig. 19

Comparative curves among experimental, BBM and MBBM for stress–strain and volume change at confining pressure 50 kPa

At a constant total confining stress, the effect of the suction is opposite to that of confining stress, since an increasing suction dilatancy. In BBM implemented in FLAC software, the trend follow the Cui and Delage [3]. experiments before the peak, however after that there are increasing dilation between two scatters. This mismatching is significant for higher suction of 1500 kPa. By plotting the results obtained from MBBM programming in FLAC, similarities to experimental results are observed. In spite of a little discrepancies in axial strains less than 5% which are probably related to the difference between laboratory and numerical analyses, but a significant number of failures in non-modified model, has been resolved.

Following their experiments and for a higher cell pressure (s₃ = 100 kPa, Fig. 20), peaks are less noticeable and appear only a little above s = 800 kPa at larger axial strains, demonstrating a lower fragility. Maximum deviator stresses increase with suction, there is only one volume change curve in this case; at s = 1500 kPa the change from contraction to dilation behavior occurs somewhat before the peak, as took place before. Comparing the results of BBM and MBBM analytical simulation for 100 kPa confining stress, proves that modified model with m factor has a good agreement with the experimental results.

Fig. 20

Comparative curves among experimental, BBM and MBBM for stress–strain and volume change at confining pressure of 100 kPa

Under s₃ = 400 kPa (Fig. 21), there are no peaks in the deviator stress curves and dilatancy disappears. The deviator stresses increase with suction whilst contraction decreases, as observed in the recent study. BBM and MBBM numerical results are in good agreement for suction s = 800 kPa whereas in higher suction this matching occurs in MBBM state. As observed in this figure, the steady state condition of stress–strain and volume change curves in all three cases are no longer available.

Fig. 21

Comparative curves among experimental, BBM and MBBM for stress–strain and volume change at confining pressure of 400 kPa

Conclusions

An elasto-plastic constitutive model based on BBM for mechanical behavior of unsaturated soils has been implemented into FLAC software. The model has been tested using a number of simulations, with regard to its implementation using the FLAC^2D user defined model capability. Also, analytical solutions were performed to calibrate the finite difference method program. The verified model shows a promising capability to study practical engineering programs, such as the response of earthen structure subjected to suction variations. Results of this study indicate that the height of structure have a remarkable effect on suction-induced settlements. A large embankment has more values of settlements during suction alternation. Therefore, in tall earthen structures potential of matric suction should be considered in design procedure. In addition to that BBM implementation in FLAC, confirm original behavior of unsaturated soils in a soil model affected by moisture changes. Based on Investigations on theorical, numerical and experimental researches, it was concluded that stress paths in LC yield curve, anisotropy and apparent tensile strength plays an effective role on mechanical and hydromechanical behavior of unsaturated natural soils which are under the influence of high suction. On the other hand, with more accurate evaluation of Alonso’s model, the following concluding remarks are drawn: (1) In Alonso’s model, in case of decreasing matric suction, LC yield curve moves that causes volume change and this makes the model to be dependent of stress path and reciprocally, when matric suction increases, the model would be stress path independent. (2) In, all the experimental studies performed by Alonso et al. [2] before saturation stage, predicted volume change is smaller than the actual one. This difference can be decreased in two ways. One is alternations of volume change coefficients, which have a significant effect on volume change after saturation that causes errors. The other is loading collapse modification in which a factor is multiplied by LC equation. This coefficient could be inversely related to suction which have a capability to control the size of plasticity zone in LC curve.

This modification improved Alonso’s model insofar as results obtained from MBBM is close to experimental results. Settlement of the embankment model with this new relation indicates better match with reality and the experimental results. For better verification of modified model, an experimental study in the literature was investigated. Comparative charts in various cell pressures confirm that modified BBM has an acceptable agreement with the experiment, especially in higher suctions which are compatible with real and natural circumstances. Moreover, the FLAC^2D, with BBM were also tested on a real model involving the interaction of multiple components of compacted soils, and embankments. Simulating these columns in unsaturated soils, appropriate results were obtained in settlement control, in different stages of depth. Results achieved in unsaturated condition without modification in BBM relations for embankment construction procedure system have been compared with MBBM state. Therefore, the Barcelona Basic Model implemented into FLAC^2D is now relatively applicable to be applied for settlement control in earthen structures for unsaturated soils. Future work will include an extension in FLAC to more complete modeling of expansive clays.