A two stage model for moisture-induced deformations in expansive soils

Abstract

Background

Moisture-induced suction changes in expansive soils due to infiltration and evaporation result in failure of civil infrastructure. The objective of this paper was to develop a two stage deformation model by simultaneously calculating soil suction and stress state. The model predictions were validated using a one-year field monitoring data.

Results

Deformations in expansive soils closely match cyclical suction changes corresponding to seasonal weather variations. Volume changes fluctuated close to the ground surface and gradually decreased with depth (overburden pressure) due to isolation from meteorological effects. The top 2 m depth was found to be the active zone susceptible to moisture variations.

Conclusions

The model agreed well with the monitoring data trends with deviations attributed to analytical assumptions in the equations, ineffective capture of antecedent soil conditions, possible soil heterogeneity and anisotropy, and hysteresis in soil saturation and desaturation.

Background

Expansive soils exhibit volumetric deformations when their water content changes. The moisture flux is governed by climatic conditions prevalent in an area. In arid and semi-arid regions of the globe, where most expansive soils occur, the moisture loss through evaporation generally exceeds the moisture gain through precipitation thereby rendering the soil unsaturated. Consequently, the soil undergoes swelling due to rainfall and shrinkage between successive rainfall events. The effect of these alternate movements exhibits differently under various types of covers. For example, longitudinal cracks in highway embankments (Jotisankasa et al. [2011]), heaved domes in residential basements (Ito and Azam [2013]), settlement in the vicinity of trees (Driscoll [1983]), and vegetation induced movements in pipelines (Clayton et al. [2010]). The cost of infrastructure repair is comparable to other natural disasters such as the estimated $1 billion in the United States of America (Phanikumar and Sharma [2006]).

Few investigations have been reported on the field behavior of expansive clays. Zhan et al. ([2007]) analyzed infiltration through a grassed expansive soil slope in Zaoyang, Hubei, China that was subjected to artificial rainfall. Likewise, Fityus et al. ([2004]) summarized a seven year monitoring data highlighting the effect of soil covers and trees on volume changes in expansive soils at a test site near New Castle, Australia. Finally, Hu et al. ([2010]) presented a one-year field data at a water distribution pipe site in the expansive soil deposit of Regina, Saskatchewan, Canada. The scarcity in reported case histories is because such undertakings are time consuming, labour intensive, and quite costly. The above-mentioned issues can be addressed through the alternative approach of numerical modeling.

Roscoe et al. ([1958]) developed the Original Cam Clay Model (OCCM) for saturated soils. Utilizing the constitutive stress–strain relationships, the elasto-plastic constitutive model consists of the following pillars: (i) elastic property; (ii) yield surface; (iii) plastic potential; and (iv) hardening rule. Roscoe and Burland ([1968]) subsequently modified the original model by incorporating elastic shear strains and an improved yield surface. To capture the behavior of unsaturated soils, Alonso et al. ([1990]) developed the Barcelona Basic Model (BBM) by correlating soil hardening with suction. Alonso et al. ([1999]) extended the above model to expansive soils by including a conceptual fabric representation: elastic deformations in microstructure and plastic deformations in macrostructure (Gens and Alonso [1992]). The Barcelona Expansive Model (BExM) was validated using laboratory test data to understand the combined effect of vertical stress and soil suction on net volumetric strains. Farulla et al. ([2007]) confirmed that the BExM model predicts the behaviour of fissured clays possessing dual porosity. These authors concluded that successive swell-shrink cycles gradually reduce plastic strains thereby converting the soil to an elastic medium.

Vu and Fredlund ([2004]) developed an elasticity based vertical displacement model for unsaturated expansive soils. The nonlinear strain was accounted for by using the coefficient of compressibility as a function of net normal stress and soil suction. The model is appropriate for predicting monotonic paths of the above two stress state variables. Adem and Vanapalli ([2013]) proposed the Modulus of Elasticity Based Method (MEBM) for deformation prediction in similar soils by using a semi-empirical estimation of the modulus of elasticity (Vanapalli and Oh [2010]) along with transient changes in soil suction (Wilson [1997]). The overestimated vertical displacement in this model is attributed to the one-dimensionality and considering soil suction as the only stress state variable.

Derived from Ito and Hu ([2011]), this study develops a two-dimensional model (by including lateral volume changes in calculating vertical deformations) that captures normal stress (overburden soil pressure) and soil suction (net atmospheric flux). The current work is based on a one-dimensional soil atmosphere model where the swelling potential was determined using laboratory test results (Azam and Ito [2012]). The two-dimensional nature of the present model helps improve the capture of actual field conditions. This was confirmed through model validation using field monitoring data, reported by Hu et al. ([2010]).

The objective of this paper was to develop a two stage model for moisture-induced deformations in expansive soils. This included a coupled soil-atmosphere model for seepage through the soil that generated suction data and suction-based displacement model that takes into account the effect of overburden soil pressure. The predicted results were validated using a one-year field instrumentation data.

Results and discussions

Figure 1 gives the variation in soil suction at selected depths with respect to time. At the ground level, the exposed surface showed median suction values of 2000 kPa (variation of ± 500 kPa) from May through August, 1500 kPa (variation of ± 500 kPa) from September through October, 1400 kPa (negligible variation) from November through February, and 1400 kPa (variation of ± 200 kPa) from March through April. These cyclical suction changes correlated well with the seasonal weather conditions, that is, multiple rainfall events during summer and fall, snow accumulation during winter, and snowmelt and rainfall during spring. Conversely, the median suction values for the covered surface were 1600 kPa (variation of ± 100 kPa) from May through August and remained at about 1400 kPa for the rest of the year. The significantly lower suction variations correspond to the isolation of the soil from seasonal weather patterns. At 2 m depth, the exposed surface showed median suction values of 1600 kPa (variation of ± 50 kPa) from May through August, 1500 kPa (variation of ± 50 kPa) from September through October, 1350 kPa (negligible variation) from November through April. Conversely, the median suction values for the covered surface were 1300 kPa (variation of ± 50 kPa) from May through November and remained at about 1250 kPa for the rest of the year. The suction decrease and the low suction variation are attributed to steady-state moisture flow (the effect of seasonal weather variations is reduced) at increased depth (Lu and Likos [2004]).

Figure 1

Suction variation versus time (modified after Ito and Hu,[2011]).

Figure 2 gives the effect of weather conditions on soil deformation at selected depths (0 m, 0.5 m, 1 m, and 2 m) during the year for exposed ground. The average vertical soil deformations corresponded well with precipitation, evaporation, and net flux data. The surface layer at an initially unsaturated state (2000 kPa ± 500 kPa) imbibed water made available from May through August by rainfall events and lost that water between consecutive rainfall events. The maximum daily heave and settlement were predicted to be 4 mm and −1 mm, respectively. The higher daily heave is attributed to the relative ease of downward seepage compared to the more energy intensive evaporation process causing settlement (van de Griend and Owe [1994]). The corresponding daily deformations remained similar during September through October because of similar variation in suction. The ground movements were negligible during winter (November through February) because precipitation was primarily received as snow that did not flow through the soil and covered the soil from any possible evaporation. Soil deformations were found to be in the order of 2 mm to −1 mm (corresponding to suction variations of 1400 kPa ± 200 kPa) during March through April because of snowmelt and rainfall events. The fluctuations in daily soil deformations gradually decreased with depth and were found to be negligible at 2 m below the ground. This progressive reduction in deformation is attributed to the overlying soil providing a cover (isolation from climatic conditions) and overburden load (cancelling the swelling movement). The model predictions indicate that the 2 m depth is the active zone where soil deformations are governed by change in atmospheric conditions.

Figure 2

Predicted soil deformations under exposed surface.

Figure 3 illustrates the influence of cyclic weather variations on soil deformation at selected depths (0 m, 0.5 m, 1 m, and 2 m) during the year for covered ground. As expected, the effect of precipitation, evaporation, and net moisture flux on average vertical soil deformations was subdued because the overlaying soil layer that covered the lower layers of the deposit reduced moisture intake. The resulting maximum daily heave and settlement were predicted to be 0.7 mm and −0.3 mm, respectively.

Figure 3

Predicted soil deformations under covered surface.

Figure 4 gives the cumulative vertical soil displacement at the surface versus time. The model predicts that the exposed soil settled by 10 mm from May through August with spikes corresponding to weather conditions. The overall settlement during the summer months is attributed to the net water deficit in the area. This trend was reversed in September through October due to several heavy rainfall events where the soil bounced back to almost the original ground level. The alternate vertical displacement during the fall months were predicted to be ±10 mm. Insignificant soil movements occurred from November through February due to the above mentioned reasons. The cyclic vertical displacements were reactivated from March through April by spring snow melt and rainfall (causing a heave of up to 3 mm) and evaporation causing shrinking of up to 5 mm. Similar soil vertical displacements albeit at a reduced level was predicted for the covered surface.

Figure 4

Comparison of modeled and monitored data.

The model predictions generally correlated well with the measured trends in data at the investigated site (obtained from Hu et al. [2010]). The variations between the predicted and the measured data are attributed to the following: (i) analytical assumptions in the equations for calculating net flux; (ii) ineffective capture of antecedent weather conditions by the model; (iii) possible soil heterogeneity, anisotropy, and vegetation at the investigated site; (iv) hysteresis in soil saturation and desaturation during the year; and (v) lack of meterological data at the investigated site.

Conclusions

Moisture-induced deformations in expansive soils are governed by the net flux at the soil-atmosphere interface and by the suction regime within the soil deposit. A two-stage two-dimensional model was developed. The model predictions generally correlated well with the trends in monitoring data at the investigated site. Deformations in expansive soils closely match cyclical suction changes corresponding to seasonal weather variations. Volume changes fluctuated widely close to the ground surface and gradually reduced with depth due to isolation from meteorological parameters and overburden pressure. The top 2 m depth was found to be the active zone susceptible to moisture variations. The variations between the predicted and the measured data are attributed to analytical assumptions in the equations for calculating net flux, ineffective capture of antecedent soil conditions by the model, possible soil heterogeneity, anisotropy, and vegetation at the investigated site, hysteresis in soil saturation and desaturation during the year, and lack of meteorological data at the investigated site. Clearly, the current model depends on an effective capture of site conditions and material properties.

Methods

Modeling process

Figure 5 describes the soil-atmosphere/soil-deformation modeling process. First, a two-dimensional soil-atmosphere model required two input categories, namely: soil properties and atmospheric parameters. The soil properties consisted of SWCC and hydraulic conductivity functions. The SWCC was estimated by using SoilVision software based on measured index properties (w, Gs, and γ _d ) and a grain size distribution (GSD) curve. Similarly, using the saturated hydraulic conductivity k _s  = 4.22 × 10⁻⁹ m/s (Shah [2011]), the estimated hydraulic conductivity function (k(ψ)) was calculated. The climate data provided the infiltration and exfiltration fluxes to determine the boundary condition and then the net flux was calculated in relation to soil properties to obtain a suction profile. The complete theoretical details of the soil-atmospheric modeling process are described by Azam and Ito ([2012]).

Figure 5

Modeling process (modified after Ito and Hu,[2011]).

The soil deformation behaviour relative to the change in atmospheric conditions was predicted by using the concept of two independent stress state variables: net normal stress, (σ - u _a ) and matric suction, (u _a - u _w ), where σ is the total stress, u _a is pore-air pressure, and u _w is pore-water pressure (Fredlund and Morgenstern [1977]). The two-dimensional soil deformation analysis required soil elasticity parameters, namely: E, soil elasticity with respect to net normal stress and H, soil elasticity with respect to matric suction (Fredlund and Rahardjo [1993]). The volume change index C _s (obtained from the net normal stress plane) was obtained from the consolidation curve (Ito and Azam [2009]). Likewise, C _m (obtained from the matric suction plane) was obtained from the SWCC and the shrinkage curve (Ito and Azam [2013]). The elasticity parameters for both the net normal stress plane and the matric suction plane were calculated using a Poisson’s ratio μ (estimated from the typical values of similar soils) for the soil structure as follows (Vu and Fredlund [2004]):

$$\mathit{E}=\frac{4.605\left(1+\mathit{\mu }\right)\left(1-2\mathit{\mu }\right)\left(1+{\mathit{e}}_0\right)}{{\mathit{C}}_{\mathit{s}}}\left(\mathit{\sigma }-{\mathit{u}}_{\mathit{a}}\right)$$

(1)

$$\mathit{H}=\frac{4.605\left(1+\mathit{\mu }\right)\left(1+{\mathit{e}}_0\right)}{{\mathit{C}}_{\mathit{m}}}\left({\mathit{u}}_{\mathit{a}}-{\mathit{u}}_{\mathit{w}}\right)$$

(2)

The net normal stress acting on a model profile before applying the weather conditions is the total vertical stress resulted from the overburden pressure. The applied atmospheric conditions adjust the void ratio along with the matric suction plane and the maximum variations in water content from full saturation to desiccation gives about 5 kPa/m stress difference, but it is almost equalized by daily fluctuations. Due to negligible variation in net normal stress state induced by the water content change, the model solely captured the simultaneous cancellation of swelling with the overburden pressure.

Soil deformations in vertical and horizontal directions was solved in accordance with the following governing equation relative to the soil suction change (Fredlund and Vu [2003]):

$$\frac{\partial }{\partial \mathit{x}}\left\{\mathit{C}\left[\left(1-\mathit{\mu }\right)\frac{\partial \mathit{u}}{\partial \mathit{x}}+\mathit{u}\frac{\partial \mathit{v}}{\partial \mathit{y}}-\frac{\left(1+\mathit{\mu }\right)}{\mathit{H}}\left({\mathit{u}}_{\mathit{a}}-{\mathit{u}}_{\mathit{w}}\right)\right]\right\}+\frac{\partial }{\partial \mathit{y}}\left\{\mathit{G}\left(\frac{\partial \mathit{v}}{\partial \mathit{x}}+\frac{\partial \mathit{u}}{\partial \mathit{y}}\right)\right\}=0$$

(3)

$$\frac{\partial }{\partial \mathit{x}}\left\{\mathit{G}\left(\frac{\partial \mathit{v}}{\partial \mathit{x}}+\frac{\partial \mathit{u}}{\partial \mathit{y}}\right)\right\}+\frac{\partial }{\partial \mathit{y}}\left\{\mathit{C}\left[\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}+\left(1-\mathit{\mu }\right)\frac{\partial }{\partial \mathit{y}}-\frac{1-\mathit{\mu }}{\mathit{H}}\left({\mathit{u}}_{\mathit{a}}-{\mathit{u}}_{\mathit{w}}\right)\right]\right\}+\mathit{pg}\mathit{=}0$$

(4)

where, p is density of soil, g is acceleration due to gravity, C = E/(1-2 μ)(1 + μ), and G = E/2(1 + μ).

The model consisted of a homogeneous soil deposit having a 6 m wide x 4 m deep geometry. To understand the effect of meteorological conditions on the soil deposit, one half of the top boundary was treated as the exposed surface (closely mimicking the vegetated park) and the other half was considered to be the covered surface (pertaining to the asphalt-pavement) at the investigated site. Zero water flux was applied at the bottom boundary to represent no ground water table at 15 m depth at the investigation site. Daily climate data from May 1, 2009 to April 30, 2010 was obtained from the Regina International Airport weather station (located about 3 km from the investigation site). An initial matric suction of 1600 kPa based on the field measurement reported by Vu et al. ([2007]) was used. In the soil-deformation model, free movements of soil in the vertical direction were allowed at the top boundary and the horizontal directions were fixed. Likewise, the lower boundary was fixed in both directions.

A general purpose partial differential equation solver, FlexPDE was used for the analysis. The solver utilized the finite element method generating a triangular mesh over a two dimensional geometry. The adequacy of the mesh was constantly calculated and the automatic mesh refinement feature of the solver was applied to reduce error within a tolerance of 0.001. Using the script editor, the governing equation and material properties were directly input in the model compiler. An equation analyzer expanded the defined equation and the material properties, performed spatial differentiation, and applied integration by parts thereby reducing second order terms to create symbolic Galerkin equations for use in the weighted residual method. The Galerkin equations were further differentiated to form the Jacobian coupling matrix for improved convergence. Likewise, the solution curvature was also calculated to include time integration for better accuracy. The model outputs were in the form of soil suction (kPa) and vertical deformation (mm) for the soil-atmosphere model and the soil-deformation model, respectively.

Material properties

Borehole drillings were made in July 2006 and May 2009 at Cross place, Regina, and the materials used for this numerical modelling had been investigated by Azam et al. ([2013]). Figure 6(a) shows the grain size distribution curves for a sample obtained in May 2009. To avoid high complexity of numerical modelling, soil properties of a representative sample were used to represent a clay deposit in the model.

Figure 6

Material properties: (a) grain size distribution, (b) soil water characteristics curve, (c) hydraulic conductivity function.

Figure 6(b) indicates the estimated SWCC. The curve exhibited the desaturation typical of fine-grained soils comprises of gradual decrease in the gravimetric moisture content as soil suction increases due to the presence of small pore spaces. Figure 6(c) shows the hydraulic conductivities as a function of soil suction. The curve also well captured the pore-size distribution derived from the grain size distribution curve shown in Figure 6(a). The reduction in hydraulic conductivity occurred slowly for the investigated clay and it only started to deviate from the saturated hydraulic conductivity value at a suction value around 100 kPa. Table 1 presents the soil properties used for the soil-deformation analysis. The elasticity parameter functions for the soil were calculated from equation [1] and [2] by using the volume change indices C _s (from net normal stress plane) and C _m (from matric suction plane) and an assumed Poisson’s ratio.

Table 1 Summary of soil properties

Climate data

Figure 7 shows the meteorological parameters used for this numerical model from May 1, 2009 to April 30, 2010. With minor variations, the data generally followed the long-term average values for the area that is classified as a semi-arid zone (BSk) according to the Köppen Climate Classification System (McKnight and Hess [2001]). Figure 7(a) presents daily precipitation along with the cumulative values. The total precipitation was recorded to be 288 mm. Rainfall primarily occurred during the summer and fall months with sporadic storm events. The accumulated snow (12 mm) gradually melted in early spring and the process was completed on March 27, 2010, when the mean temperature peaked to 1.4°C and remained above freezing for most of the day. Figure 7(b) shows the daily mean temperature. The annual mean temperature was found to be 2.7°C. The average from May to September was 14.4°C and the average from October to April was −5.8°C. The maximum difference between the high and low temperatures was 57.9°C. Relative humidity and wind speed data are shown in Figure 7(c) and (d), respectively. The relative humidity data was generally inversely proportional to the temperature data. The average of relative humidity in winter was 78% and reduced to 64% in the summer. There was no clear trend for wind speed – it changed regardless of the season. The daily net radiation was estimated in Figure 7(e) from the solar radiation measurement from November 1, 1999 to October 31, 2000 recorded in Swift Current, SK, by using albedo values of 0.1 for the period from November 1 to March 31 and 0.3 for the rest of the year.

Figure 7

Climate data: (a) precipitation; (b) temperature; (c) relative humidity; (d) wind speed; (e) net radiation.

The above data implies that during the summer time, the average temperature remained high (14.4°C), the average wind speed was consistent (17.8 km/hour), the average relative humidity was low (64%), and the average net radiation was high (13.9 MJ/m²/day). Such conditions generally render the surface layer of Regina clay desiccated. Therefore, in any sporadic rainfall event during the summer, the expansive clay can readily adsorb all of the available water thereby resulting in swelling.