3 Coupled Consolidation and

Solute Transport Model:

1D Finite Deformation
3.1 INTRODUCTION
Land-based containment facilities are commonly used for the disposal of municipal
solid waste and contaminated dredged material (Liu, 2007). In modern landfills,
liner systems are designed to isolate the landfill contents from the surrounding
environment to protect the groundwater from pollution. For well-constructed
composite liners, the geo-membrane typically has few defects, so restricting
advection through it (Giroud and Bonaparte, 1989; Foose et al., 2002). However,
volatile organic compounds (VOCs) can diffuse through membranes with magnitude
four to six orders greater than the possible advection. Therefore, diffusion of VOCs
in composite liners is viewed as a critical issue in the design of landfill liners (Foose,
2002).
The VOC transit time was traditionally estimated using the diffusion equation
(Rowe and Badv, 1996; Fityus et al., 1999; Foose, 2002). However, several field
tests have reported that the transit of VOCs is much earlier than theoretical
predictions (Workman, 1993; Othman et al., 1997). Many researchers attribute
this to consolidation and associated advective transport. Several theoretical models
coupling mechanical consolidation with solute transport were constructed in recent
years (Smith, 2000; Fox, 2007b; Lewis et al., 2009).
Based on the 1D Biot consolidation theory, Zhang et al. (2012a) proposed
an advection-diffusion equation that incorporates the degree of saturation,
compressibility of the pore fluid (CPW) and dispersivity of the solute transport in a
nearly saturated deforming porous medium. Both CPW and dispersivity were found
to influence solute migration within the CCL, significantly so in some circumstances.
However, Zhang et al. (2012a) considered an infinitesimal strain, (i.e., small
deformation) model. Additionally, they did not consider the material and geometric
nonlinearity, factors that could be important in some circumstances (Lewis et al.,
2009). Financial constraints sometimes limit deployment of the relatively costly
CCLs. Natural clay deposits (sometimes with relatively high compressibility) are
used as substitutes. Since the soft clayey soil generally provides a good contact
adhesion with a geomemebrane, high effectiveness is a priori expected. However,
the finite deformation caused by the emplacement of waste cannot be neglected.
In this chapter, the small deformation model for solute transport in a nearly
saturated medium (Zhang et al., 2012a) will be further extended to finite
deformations. This allows us to clarify the influence of consolidation in the progress

49
50 Poro-Elastic Theory with Applications to Transport in Porous Media

of solute transport (using a time-dependent boundary in terms of void ratio at the

CCL base). The influence of the degree of saturation on the VOC transit time in clay
barriers will also be examined. To account for the geometric non-linearity, a material
coordinate system is used. Both CPW and dispersivity are considered in the new
model. Further, our approach incorporates nonlinearity of the constitutive properties
related to soil compressibility, the hydraulic conductivity and decreasing effective
diffusion coefficient. A parametric study is carried out to examine the influence of
several dominant parameters on the process of solute transport in porous medium.

3.2 MODEL FORMULATION

Lewis et al. (2009) and Peters and Smith (2002) developed a model coupling finite
strain consolidation and solute transport in a fully saturated soil. Below, the CPW
and dispersion in a nearly saturated soil is included.

3.2.1 FINITE STRAIN CONSOLIDATION

A Lagrangian coordinate system (z, t) is employed to derive the flow and transport
equations. We define ξ (z,t) as the particle displacement with ξ (z, 0) = z. The
relationship between Lagrangian and Eulerian (ξ ,t) coordinate systems then implies
that for any variable F(z,t) = f (ξ (z,t),t):

∂F ∂ f ∂ξ ∂F ∂ f ∂ξ ∂ f ∂f ∂f
= , = + = vs + , (3.1)
∂z ∂ ξ ∂ z ∂t ∂ ξ ∂t ∂t ∂ξ ∂t

where vs = ∂ ξ /∂t is the solid velocity.

The equation describing changes in void ratio, e(z,t), are derived from the
continuity equations for the solid and fluid phases together with Darcy’s law. The
mass balance equation of the solid phase in differential form is
 
∂ ∂ξ
ρs (1 − n) = 0, (3.2)
∂t ∂z

where ρs is the soil grain density, n = e/(1 + e) is the current porosity, and n0 =
n(z, 0) is the initial porosity. Note that, for constant ρs , the Jacobian, M, for the
coordinate transformation is
∂ξ 1 − n0 1+e
M= = = , (3.3)
∂z 1−n 1 + e0
where e0 is the initial void ratio.
The continuity equation for the fluid phase (i.e., pore-water) is
 
∂ ∂ξ ∂
nSr ρ f = − (ρ f q), (3.4)
∂t ∂z ∂z

where ρ f is the pore fluid density.

1D finite strain coupled model for consolidation and solute transport 51

According to Darcy’s Law, the fluid flux is given by

kv ∂ p
q=− , (3.5)
ρf g ∂ξ
where kv is hydraulic conductivity and p is excess pore pressure. If the hydraulic
gradient is constant, the Darcy equation in terms of total pressure can be transformed
to this form (Peters and Smith, 2002).
Assuming ρ f varies with pore pressure as ∂ ρ f /∂ p = β ρ f (Barry et al., 2007),
substituting (3.5) into (3.4), then the continuity equation for the fluid phase becomes:

∂ξ ∂ p ∂ 1 ∂ ∂ p ∂z
   
∂ξ
nSr β + Sr = kv , (3.6)
∂ z ∂t ∂t ∂z ρf g ∂z ∂z ∂ξ
where the compressibility of pore fluid (β ) is defined in (2.2).
Because n and n0 (implicitly embedded in ∂ ξ /∂ z) appear simultaneously, and n
is unknown, (3.6) can not be directly solved in terms of p. In the following derivation,
it turns out that once the relationship between the derivative of p (with respect to t
and a) and the corresponding derivative of e is known, it is straightforward to convert
(3.6) to an equation in terms of e.
Assuming self-weight is negligible due to the relatively small thickness of the
CCL (Zhang et al., 2012a), the vertical force equilibrium is
∂σ
= 0, (3.7)
∂z
where σ (now a function of t only) is the total normal stress of the soil and the z
coordinate is vertically upwards. Assuming the compressive normal stress is positive,
i.e., σ = σ 0 + p (σ 0 is the effective normal stress), (3.7) leads to:
∂p ∂
∂z 1 + e0 1 ∂ e
= −σ 0 + σ = , (3.8)
∂ξ ∂z ∂ξ 1 + e αv ∂ z
where αv = −de/dσ 0 is the coefficient of soil compressibility.
In the absence of self-weight, the rate of change of total stress at an arbitrary
location equals that of the external top loading,
∂σ ∂Q
= , (3.9)
∂t ∂t
where Q is the external load. The rate of change of the excess pore water pressure in
the time domain is
∂p ∂ ∂Q 1 ∂e
= (σ − σ 0 ) = + . (3.10)
∂t ∂t ∂t αv ∂t
Substituting (3.3), (3.8), (3.10) into (3.6) yields:
eSr β Sr ∂ e 1 + e0 ∂ kv ∂e Sr β e ∂ Q
   
+ − =− . (3.11)
(1 + e0 )αv 1 + e0 ∂t ρ f g ∂ z αv (1 + e) ∂ z 1 + e0 ∂t
52 Poro-Elastic Theory with Applications to Transport in Porous Media

For the fully saturated case and when the CPW is neglected, i.e., β = 0, (3.11)
reduces to:
1 ∂ e 1 + e0 ∂ kv ∂e
 
= , (3.12)
1 + e0 ∂t ρ f g ∂ z αv (1 + e) ∂ z

which is identical to Eq. (1) of Lewis et al. (2009).

3.2.2 SOLUTE TRANSPORT EQUATIONS

Solute transport occurs in both solid and fluid phases. Here, for the nearly saturated
soil, the mixture of pore-water and entrapped air is taken as a homogeneous fluid.
Due to the discrete air bubbles, VOC transport by gas diffusion can be neglected in a
nearly-saturated soil. Therefore, the mass conservation equation for the solute in the
solid phase is
 
∂ ∂ξ 0
(1 − n) ρs S = fa→s , (3.13)
∂t ∂z

where S is the mass of solute sorbed on or within the solid phase per unit mass of
the solid phase and fa→s
0 denotes rate of solute loss in the water phase by solid phase
sorption.
The mass conservation equation for solute in the fluid phase is

∂ Jf
 
∂ ∂ξ 0
nSr c f =− − fa→s , (3.14)
∂t ∂z ∂z

where c f is the concentration of the solute in the pore fluid. In (3.14), the term ∂ ξ /∂ z
comes from the volumetric change (Peters and Smith, 2002) and J f represents solute
flux in the fluid phase, which is described by (Peters and Smith, 2002):

nSr D ∂ c f
J f (z,t) = nSr (v f − vs )c f − , (3.15)
M ∂z
where D is the hydrodynamic dispersion coefficient. It is given by the sum of the
effective diffusion coefficient (De ) and the coefficient of mechanical dispersion (Dm ):

Dm = αL v f − vs , (3.16)

where αL is dispersion coefficient, v f is the pore fluid velocity and v f − vs denotes

the relative velocity of the pore fluid.
Based on (3.13)–(3.15), we have:

∂ nSr D ∂ c f
   
∂   ∂ξ ∂ 
nSr c f + (1 − n)ρs S nSr (v f − vs )c f . (3.17)

= −
∂t ∂z ∂z M ∂z ∂z
1D finite strain coupled model for consolidation and solute transport 53

The above equation can be further simplified with Darcy’s Law, (3.5), and the
mass balance equations for both solid and fluid phases, (3.2) and (3.4), respectively.
(3.17) can then be expressed as:

∂ξ ∂cf ∂ ξ ∂ S ∂ nSr D ∂ c f kv ∂ p ∂ c f
 
nSr + (1 − n)ρs = +
∂ z ∂t ∂ z ∂t ∂ z M ∂z ρf g ∂ξ ∂z
(3.18)
∂ ξ ∂ p β kv ∂ p ∂ p
 
+ nSr β − cf .
∂ z ∂t ρ f g ∂ ξ ∂ z

Substituting (3.8) and (3.10) into (3.18) results in:

e ρs Kd ∂ c f ∂ e(1 + e0 ) ∂ c f
   
Sr + = Sr D
1 + e0 1 + e0 ∂t ∂ z (1 + e)2 ∂ z
kv 1 + e0 ∂ e ∂ c f
+
ρ f g αv (1 + e) ∂ z ∂ z

e

∂Q 1 ∂e
 (3.19)
+ β Sr +
1 + e0 ∂t αv ∂t
 2 #
kv 1 + e0 ∂ e
− cf ,
ρ f gαv2 1 + e ∂ z

where Kd describes the partitioning coefficient.

3.2.3 SPECIAL CASES

In this section, three special cases of the present model are outlined.

A. Saturated soil with finite deformation

For a saturated soil, where Sr = 1, and incompressible pore fluid, i.e., β = 0, (3.19)
reduces to:
e ρs Kd ∂ c f ∂ e(1 + e0 ) ∂ c f kv 1 + e0 ∂ e ∂ c f
   
+ = 2
D + , (3.20)
1 + e0 1 + e0 ∂t ∂ z (1 + e) ∂z ρ f g αv (1 + e) ∂ z ∂ z

which is identical to Eq. (4) of Lewis et al. (2009) and Eq. (44) in Peters and Smith
(2002).

B. Small deformation model

Under the assumptions of negligible self-weight and small deformation (constant
porosity, i.e., n = n0 ), the coupled deformation model can be expressed the following,
which is identical to that of Zhang et al. (2012a):

∂p ∂ 2u 1 ∂ ∂p
 
Sr n0 β + Sr = kv , (3.21a)
∂t ∂t∂ ξ ρw g ∂ ξ ∂ξ
54 Poro-Elastic Theory with Applications to Transport in Porous Media

2(1 − ν) ∂ 2 u ∂p
G = , (3.21b)
(1 − 2ν) ∂ ξ 2 ∂ξ
and
∂cf ∂ 2c f kv ∂ p ∂ 2 c f
[Sr n0 + (1 − n0 ) ρs Kd ] = Sr n0 De 2
− αL
∂t ∂ξ ρw g ∂ ξ ∂ ξ 2
∂cf ∂p ∂ 2u

+ −αL Sr n0 β − αL Sr
∂ξ ∂t ∂ ξ ∂t
 2
αL β kv ∂ p ∂ 2u
+ + Sr De (1 − n0 ) 2
ρw g ∂ξ ∂ξ
(3.21c)
kv ∂ p ∂u

+ − [Sr n0 + (1 − n0 ) ρs Kd ]
ρw g ∂ ξ ∂t
 2
∂p kv ∂p
+ Sr n0 β cf −β cf
∂t ρw g ∂ ξ
∂u ∂ p
+ Sr n0 β cf ,
∂t ∂ ξ
where u is the soil displacement, G is the shear modulus and ν is Poisson’s ratio.
The constant material coefficients can be described as:
cv ρ f g(1 − 2ν) (1 + e p )(1 − 2ν)
G= = ,
2kv (1 − ν) 2(1 − ν)αvp (3.22)
kv = k p , De = De0 ,

where cv is the consolidation coefficient; ks and k p the saturated hydraulic

conductivity and hydraulic conductivity of the soil corresponding to e p (the void
ratio corresponding to pre-consolidation stress), respectively.

C. Nearly saturated soil with no deformation

For the partially saturated no deformation model, i.e., e = e0 , ξ = z, the overloading,
Q, does not affect solute transport. In the spatial coordinate system (ξ , t), (3.19)
reduces to the linear diffusion equation:

∂cf ρs Kd −1 ∂ 2 c f
 
= D 1+ . (3.23)
∂t Sr e0 ∂ξ2

3.3 VARIATIONS OF PARAMETERS IN CONSOLIDATION AND

SOLUTE TRANSPORT PROCESSES
The finite deformation model allows consideration of the effects of variations in the
coefficients of consolidation and transport (such as the coefficient of compressibility,
1D finite strain coupled model for consolidation and solute transport 55

αv , hydraulic conductivity, kv and hydrodynamic dispersion, D) on solute transport

process. Lewis et al. (2009) utilized void ratio-dependent functions for the related
coefficients while Li and Liu (2006) used a fractal pore-space theory to develop
fractal models of water flow and solute diffusion in rigid unsaturated soils. Their
approach allowed comparison of these coefficients between the fully saturated and
unsaturated cases. Here, a combination of both models is employed so that the
hydraulic conductivity and the effective diffusion depend on both the void ratio and
the degree of saturation. Linear, reversible solute sorption is assumed in this study.
however, the approach can be adapted for other sorption models.

3.3.1 SOIL COMPRESSIBILITY

The soil layer is assumed to be over-consolidated, and compression of the soil
layer commences when the applied stress exceeds the pre-consolidation stress, i.e.,
deformation due to re-compression is neglected. In this case, the void ratio is
idealized as a linear function of the logarithm of the effective stress (Means and
JV., 1964):
0
!
σ
e = e p −Cc log , (3.24)
σ p0
where σ 0 is effective stress, σ p0 denotes the pre-consolidation stress and Cc is the
compression index of the soil (defined by the absolute value of the slope of the
idealized virgin compression line).
For a nearly saturated soil, the degree of saturation is sufficiently high so that
the air phase exists in the form of occluded bubbles. Vaughan (2003) claimed that
the presence of occluded air bubbles is unlikely to affect soil effective stresses.
Therefore, (3.24) is employed to describe the volumetric change of a nearly saturated
soil.
The coefficient of compressibility in terms of void ratio can be obtained by
differentiation of (3.24) with respect to effective normal stress (Lewis et al., 2009):
e − ep
  
αv = αvp exp ln10 , (3.25)
Cc
0
where αvp is the coefficient of compressibility corresponding to σ p , i.e.,
Cc
αvp = . (3.26)
σ p0 ln10

3.3.2 HYDRAULIC CHARACTERISTICS

For hydraulic conductivity, an empirical relationship describing its variation with
void ratio in saturated clay soils is given as (Mitchelll, 1993):
e − ep
  
ks = k p exp ln(10) , (3.27)
Ck
where Ck is the hydraulic conductivity index.
56 Poro-Elastic Theory with Applications to Transport in Porous Media

The power law relationship equation for hydraulic conductivity versus water
content θ (= Sr n) is (Li and Liu, 2006):
 α
θ
kv = ks , (3.28)
θs
where θs is saturated water content, and α falls in the range of 2.68 to 2.78 for clay
loam.

3.3.3 DISPERSION COEFFICIENT

In a saturated soil, the effective solute diffusion coefficient is defined as the product
of the free diffusion coefficient of the solute in the pore fluid (D f ) and the tortuosity
factor (t f ), which accounts for the irregular path that diffusing molecules must take
through the pore space (Acar and Haider, 1990). Lewis et al. (2009) claimed that it
is rational to take De as constant, because uncertainty of the range of τ f can be the
same order of consolidation-induced change of De . Alternatively, the reduction of
De can be expressed with a hypothetical relationship associated with the overall void
ratio change as (Lewis et al., 2009; Morel-Seytour et al., 1996):
e0 − e e−ef
 
De = + De0 , (3.29)
3(e0 − e f ) e0 − e f
where e f denotes the final void ratio, and De0 is the initial effective dispersion
coefficient.
In variably saturated soils, the effective diffusion coefficient, De , depends on soil
water content, bulk density, and soil type for soils with different textures. Regarding
the water content, there is a threshold value under which solute diffusivity vanishes
(Hunt and Ewing, 2003; Hamamoto et al., 2009). The impedance factor (Porter
et al., 1960) (i.e., the ratio of solute diffusion coefficient in soil to product of solute
diffusion coefficient in free water and volumetric soil water content), decreased with
increasing bulk density for each soil type, but the effect of the overall bulk density
on the impedance factor is minor compared with the effect of soil water content and
soil type (Hamamoto et al., 2009).
The effective diffusion coefficient was found to decrease with decreasing
saturation in laboratory experiments (Barbour et al., 1996). The decrease was found
to be quite rapid initially, followed by a near-linear decline for degree of saturation
below 60%. Here, the soil diffusion coefficient is expressed as (Li and Liu, 2006):
De = 1.1D f θ (θ − θt ), (3.30)
where θt denotes threshold water content, which was observed to become higher
with increasing clay content and varies between 3% and 20% for clay soil.

3.3.4 SORPTION
It has been reported that the effect of the degree of saturation on the adsorption
coefficient is insignificant from full saturation to a degree of saturation of 10%
1D finite strain coupled model for consolidation and solute transport 57

(Barbour et al., 1996). A significant decrease in the adsorption coefficient only

occurs in cases with a low degree of saturation. In this study the degree of saturation
varies from 1 to 0.8, i.e., the effect on sorption can be neglected. Therefore, the
concentration of solute in the solid phase, S, is expressed as:

S = Kd c f . (3.31)

This assumption of a linear sorption is valid at the relatively low concentrations that
are usually found in the municipal waste disposal sites (Mathur and Jayawardena,
2008).

3.4 APPLICATION TO A LANDFILL LINER

As the schematic in Figure 2.1 shows, the composite landfill liner beneath a primary
leachate collect system (PLCS) consists of an impermeable (to diffusion of inorganic
solute) geomembrane, an underlying engineered compacted clay layer (CCL), and a
second leachate collecting system (SLCS).
The model parameters employed in the following analyses are based on those
used in recent studies of solute transport in composite liners (Foose, 2002; Lewis
et al., 2009). Because of the unavailability of consolidation data in the literature,
hypothetical values of the applied stress, pre-consolidation stress, compression
index, hydraulic conductivity index, threshold moisture content and other parameters
in calculating the De and kv are used. As a primary parameter, the compression index
covers a large range to account for the high-compressibility soil considered (Lewis
et al., 2009). However, the related applied stress was selected to avoid negative and
unrealistically low void ratios. The parameters used are given in Table 3.1.

3.4.1 BOUNDARY CONDITIONS FOR CONSOLIDATION

The following boundary conditions are introduced. Assuming there are no defects
in the geomembrane, the top boundary (z = 0) is assumed to be impermeable, i.e.,
q = 0. Therefore, from (3.5) and (3.8),

∂e
= 0 at z = 0. (3.32)
∂z
At the bottom drainage boundary (z = L), the excess pore pressure is zero and a
Dirichlet-type boundary condition for void ratio (e) can be derived from the effective
stress–void ratio equilibrium relationship, (3.24):
!
σL0
e = e p −Cc log , (3.33)
σ p0

where σL0 denotes the effective stress at bottom.

The excess pore pressure vanishes at the bottom boundary, so σL0 = σa , where σa
is a time-varying stress due to the external overburden. Note that σa is the maximum
58 Poro-Elastic Theory with Applications to Transport in Porous Media

Table 3.1
Values of input parameters.

Parameter Value
Maximum applied stress (ramp loading for 2 years), σa 450 kPa
Preconsolidation stress, σ p0 50 kPa
Compression index, Cc 0.2, 0.8
Preconsolidation hydraulic conductivity, k p 10−9 , 2×10−10 m/s
Constant, α 2.7
Hydraulic conductivity index, Ck 0.585
Thickness of geomembrane, h 0.0015 m
Thickness of CCL, L 1.22 m
Mass transfer coefficient of geomembrane, PG 4 × 10−11 m2 /s
Initial effective diffusion coefficient, De0 2 ×10−10 m2 /s
Free diffusion coefficient in the pore fluid, D f 10−9 m2 /s
Threshold moisture content, θt 0.05
Partitioning coefficient, Kd 0, 0.2, 1 ml/g
Dispersion, αL 0, 0.1 m
Initial void ratio, e0 ( = e p ) 1.17
Acceleration due to gravity, g 9.81 m/s2
Initial density of pore-water, ρ f 103 kg/m3
Density of the solid phase, ρs 2.7 × 103 kg/m3
Degree of saturation of clay, Sr 1, 0.9, 0.8

loading in the model of Lewis et al. (2009). The void ratio rapidly approaches a
steady value, which consequently leads to a spurious higher fluid velocity and faster
solute transport. To distinguish the cases, we label the present boundary condition at
the CCL bottom as ‘BCC’ and ‘BCL’, i.e., the boundary conditions used by Lewis
et al. (2009).

3.4.2 BOUNDARY CONDITIONS FOR SOLUTE TRANSPORT

At the top of the CCL, VOC diffusion through the geo-membrane is described by
Fick’s law (Booker et al., 1997), so the concentration gradient is proportional to the
difference in concentrations on each side of the (sufficiently thin) geomembrane. In
the material coordinate system, the boundary condition is (Lewis et al., 2009):

∂cf (1 + e(0,t))2 PG
c f (0,t) −C f 0 , (3.34)


(0,t) =
∂z e0 (1 + e0 ) hDe

where C f 0 is the (constant) solute concentration at the top surface of the

geo-membrane with the assumption that the landfill waste volume is large (Peters and
1D finite strain coupled model for consolidation and solute transport 59

Smith, 2002); h and PG are, respectively, the thickness and the permeation coefficient
for the solute in the geo-membrane.
The lower boundary condition for the solute concentration (c f ) is (Peters and
Smith, 2001):
∂cf
= 0, at z = L, (3.35)
∂z
which assumes negligible diffusion below the CCL base (Barry and Sposito, 1988).

3.4.3 MODEL VERIFICATION

Since there are no experimental data available in the literature, the present model
was reduced to the full-saturation case using the same boundary condition at the
CCL bottom for e as used by Lewis et al. (2009), i.e., σa is taken as the maximum
loading; and Kd = 0, αL = 0, Cc = 0.8, k p = 10−9 m/s. A comparison between the
present and previous models is illustrated in Figure 3.1. In the figure, the results
of the finite deformation with constant and decreasing hydrodynamic dispersion,
3.29, small deformation model (Zhang et al., 2012a) and the pure diffusion model
(i.e., no deformation model) are included. Both consolidation (i.e., void ratio, e,
distribution) and relative concentration obtained from the present model are in
excellent agreement with results of Lewis et al. (2009). As shown in Figure 3.1,
with the constant effective diffusion coefficient, the small deformation model (Zhang
et al., 2012a) predicts a slower solute migration than the corresponding finite
deformation model.

−1
10
Relative concentration at CCL base, Cf/Cf0

−2
10

−3
10 LD, Constant De
LD, Decreasing De
ND
SD
Lewis et al. (2009)
−4
10
5 10 15 20 25 30 35 40
Time, t (y)

(a) Void ratio evolution (b) Breakthrough curves

Figure 3.1 Comparison of (a) void ratio evolution and (b) breakthrough curves between the
present model (solid line) and Lewis et al. (2009) (circle). Notations: FD: finite deformation
model, SD: small deformation model, ND: no deformation model.

3.4.4 CORRECTNESS OF THE BOUNDARY CONDITION AT CCL BASE

The differences due to the different boundary conditions, ‘BCL’ (used by Lewis et al.
(2009)) and ‘BCC’ (used in the present model), are presented in Figure 3.2, where
60 Poro-Elastic Theory with Applications to Transport in Porous Media

Cc = 0.8 and k p = 10−9 m/s. A comparison of Figure 3.2(a) (BCC) and 3.1(a) (BCL)
shows that taking σa as the maximum loading leads to a greater void ratio gradient
and a faster consolidation process, although the final value of e is very close. This
initially speeds up the solute transit slightly, and then slows it down in the long-term
(Figure 3.2(b)). The reason the trend reverses after the consolidation completes for
the ‘BCL’ case is that the higher solute concentration level during the consolidation
phase of ‘BCC’ occurs later resulting in an increased advective flux. The separation is
more obvious for the relatively soft and higher permeablility cases. In the following
sections all numerical results are based on the boundary condition ‘BCC’.
0
10

Relative concentration at CCL base, Cf/Cf0

−1
10

−2
10

Case 1(BCC)
−3
10 Case 2(BCC)
Case 3(BCC)
Case 1(BCL)
Case 2(BCL)
Case 3(BCL)
−4
10
0 5 10 15 20 25 30 35 40
Time, t (y)

(a) Void ratio evolution (BCC only) (b) Breakthrough curves (BCC and BCL)

Figure 3.2 Influence of Boundary condition of void ratio (e) at CCL base (a) void ratio
evolution (BCC only) and (b) breakthrough curves (Sr = 1, β = 0, αL = 0, constant De ). In
(b), solid line for ‘BCC’, and dash-dot line for ‘BCL’. Case 1: k p = 2 × 10−10 m/s, Cc =0.8;
Case 2: k p = 10−9 m/s, Cc =0.8; and Case 3: k p = 10−9 m/s, Cc =0.2.

3.5 NUMERICAL RESULTS AND DISCUSSIONS

A numerical solution was constructed using COMSOL 3.5a (COMSOL, 2010). It
discretized the domain into unstructured Lagrange-linear elements with a maximum
global element size of 10−2 m, and maximum local element size at the end
boundaries (where the most rapid changes occur) of 10−4 m. Temporally, the
sub-time step was 10−2 y. To be easily interpreted, solution curves were plotted in
the spatial coordinate x:
Z L
e0 − e(ζ )
x = z+ dζ . (3.36)
z 1 + e0
Thus, the first-order PDE,
∂x e0 − e(z)
= 1− , (3.37)
∂z 1 + e0
1D finite strain coupled model for consolidation and solute transport 61

with boundary conditions x(0,t) = Smt and x(L,t) = L was constructed to find x,
where the settlement Smt is given by:
Z L
e0 − e(ζ )
Smt = dζ . (3.38)
0 1 + e0

3.5.1 EFFECT OF CONSOLIDATION

On basis of the ‘BCL’ boundary condition, Lewis et al. (2009) observed that there
is no noticeable solute concentration at the CCL base when consolidation of the
liner is completed even for the case of very high compressibility (Cc = 0.8). They
thus concluded that transport can be simulated using the pure diffusion model with
the final void ratio value. However, during consolidation the distribution of solute
concentration changes, which is the initial condition of what follows. Thus, advective
transport due to consolidation may not be negligible.
0 0

−0.1 −0.1
Settlement at CCL top, S (m)

Settlement at CCL top, S (m)

−0.2 −0.2

−0.3 −0.3

−0.4 −0.4

−0.5 −0.5
kp = 2×10−10 m/s, Cc=0.8 kp = 2×10−10 m/s, Cc=0.8
kp = 10−9 m/s, Cc=0.2 kp = 10−9 m/s, Cc=0.2

0 5 10 15 20 25 30 35 40 0 20 40 60 80 100 120
Time, t (y) Time, t (y)

Figure 3.3 Consolidation settlements in a (a) saturated soil (Sr = 1) and (b) partially
saturated soils (Sr = 0.8).

Figures 3.3(a) and 3.4(a) illustrate the consolidation processes and solute
transport in a saturated soil for two cases with different compression indices (Cc )
and hydraulic conductivities (kv ). Consolidation lasts 2.2 and 34.5 y for Cc = 0.2 and
Cc = 0.8, respectively. For the ‘soft’ case, a noticeable concentration difference from
the no deformation model appears at the CCL base during consolidation, as shown in
Figure 3.4(a). The difference decreases with higher levels of sorption (Figure 3.4(a)).
The effect of consolidation on transport exists during both the consolidation and
post-consolidation stages, which is consistent with Fox (2007b). Since the advection
results in a notable concentration level at the CCL base, simplifying assumptions
such as instant deformation, pure diffusion and finite deformation without advection
modeling are not appropriate. The magnitude of solute concentration C f in Figure
3.4(a) is an order greater than that in Figure 3.4(a) . Here, the influence of sorption
is noticeable as it drastically retards the solute transport.
Figures 3.3(b) and 3.4(b) present the results for a nearly saturated soil. We see
again that soft clay consolidation has a noticeable effect on solute transport (Figure
62 Poro-Elastic Theory with Applications to Transport in Porous Media

Relative concentration at CCL base, Cf/Cf0

0.6
kd = 0.2 ml/g
Relative concentration at CCL base, Cf/Cf0

−1
0.5 10

0.4
−2
10
0.3 k = 1 ml/g
d

0.2 −3
10

0.1

−4
10
0 0 20 40 60 80 100 120
0 5 10 15 20 25 30 35 40
Time, t (y) Time, t (y)

Kd = 0 Kd 6= 0
(a) saturated soil (Sr =1)
Relative concentration at CCL base, Cf/Cf0
k = 0.2 ml/g
d
−2
10
Relative concentration at CCL base, Cf/Cf0

0.5

0.4 −4
10

0.3
−6
kd = 1 ml/g
10
0.2
−8
10
0.1

−10
10
0 0 20 40 60 80 100 120
0 20 40 60 80 100 120
Time, t (y) Time, t (y)

Kd = 0 Kd 6= 0
(b) partially saturated soils (Sr = 0.8)

Figure 3.4 Effect of consolidation on relative concentration C f /C f 0 in a (a) saturated soil

(Sr = 1, , αL = 0, constant De ) and (b) partially saturated soils (Sr = 0.8, with CPW, αL = 0.1
m, varying De as in (3.30)) with two cases (Kd = 0 and (b) Kd 6= 0) . Notations: solid line (FD,
finite deformation model): Cc = 0.8, k p = 2 × 10−10 m/s; dash-dot line (FD, finite deformation
model): Cc = 0.2, k p = 10−9 m/s; and dashed line: no deformation model (ND).

3.3(b)). However, since the effective diffusion (De ) reduces with deformation,
concentrations for the pure diffusion model surpass those of coupled models, as is
obvious for the case of Kd = 1 ml/g.
Consolidation effects are composed of the variation of void ratio and the
occurrence of pore-water flow, which in turn causes the advective transport flux.
As mentioned previously, Lewis et al. (2009) claimed the advection component can
be ignored as long as the variation of void ratio is considered. Here, we included
in Figure 3.5 the case of finite deformation without advection, i.e., advection is
removed from (3.19). Exclusion of advection underestimates the concentration level
and consequently leads to a longer transit time. In the absence of sorption, at the
nominal 10% breakthrough, a nearly twofold change occurs in the transit time; this
change increases when sorption is included.
1D finite strain coupled model for consolidation and solute transport 63

Relative concentration at CCL base, Cf/Cf0 −2

Kd = 0 ml/g
10

Kd = 1 ml/g
−4
10

−6
10
0 20 40 60 80 100
Time, t (y)

Figure 3.5 Effect of advection flux on concentration level at CCL base for partially saturated
cases (Sr = 0.8, with CPW, αL = 0.1 m, varying De as in (3.30)). For finite deformation model,
solid line: Cc = 0.8, k p = 2 × 10−10 m/s; dash-dot line: without advection flux in transport,
(3.19); dashed line: No deformation model.

3.5.2 EFFECT OF DEGREE OF SATURATION

0
10
Relative concentration at CCL base, Cf/Cf0

−1
10

−2
10
Kd = 0, Sr = 0.8
Kd = 0, Sr = 0.9

−3 Kd = 0, Sr = 1.0
10 Kd = 1 ml/g, Sr = 0.8
Kd = 1 ml/g, Sr = 0.9
Kd = 1 ml/g, Sr = 1.0
−4
10
0 20 40 60 80 100 120
Time, τ (y)

Figure 3.6 Effect of saturation Sr on transport for no-deformation model.

Figure 3.6 demonstrates that the higher saturation of the no-deformation (ND)
model results in faster solute transport due to the saturation (Sr )-dependent effective
diffusion; the gap is larger in the presence of sorption. Concentrations predicted
by the coupled finite deformation and solute transport model are shown in Figure
3.7. For cases with parameters Cc = 0.8 and k p = 10−9 m/s, consolidation lasts for
approximately 12.8 y. Higher saturation results in faster solute transport because of
64 Poro-Elastic Theory with Applications to Transport in Porous Media

greater effective diffusion, regardless of the sorption. For decreasing De , the transit
time increases, as shown in Figure 3.7(a). With sorption, finite deformation with
Sr = 0.8 and constant De leads to almost the same concentration as for the ND
model (Figure 3.7(b)). Again, this demonstrates that the effect of unsaturation is
more apparent in the presence of sorption. Interestingly, with both sorption and
decreasing De taken into account, finite deformation (FD) models will not always
produce faster solute transport (Figure 3.7(a)). During consolidation and in the early
post-consolidation stage, the FD models have a faster transit, but then are surpassed
by the ND model because the effective diffusion is reduced due to compaction.
However, the decreasing De with compaction is inevitable. In the field, VOC has
been shown to appear earlier than predicted by the pure diffusion model has been
observed (Peters and Smith, 2002). Possible explanations are: (1) the constitutive
relationships for soil parameters are not accurate enough; or (2) other factors, such
as heat transfer, should be also included in the model. Relative concentration at CCL base, Cf/Cf0
Relative concentration at CCL base, Cf/Cf0

0.4
FD, Sr = 1.0
1e−2
FD, S = 0.9
r
0.3 FD, Sr = 0.8
ND, Sr= 1.0 1e−4

0.2
1e−6
FD, Sr = 1.0
0.1 FD, Sr = 0.9
1e−8
FD, Sr = 0.8
ND, Sr= 1.0
0 0
0 10 20 30 40 0 20 40 60 80 100 120
Time, τ (y) Time, τ (y)

Kd = 0 Kd = 1 ml/g
(a) The case with decreasing De
Relative concentration at CCL base, Cf/Cf0

Relative concentration at CCL base, Cf/Cf0

0.7 FD, Sr = 1.0 FD, Sr = 1.0

0.25
FD, Sr = 0.9 FD, Sr = 0.9
0.6
FD, Sr = 0.8 FD, Sr = 0.8
0.2
0.5 ND, Sr= 1.0 ND, Sr= 1.0

0.4 0.15

0.3
0.1
0.2
0.05
0.1

0 0
0 10 20 30 40 0 20 40 60 80 100 120
Time, τ (y) Time, τ (y)

Kd = 0 Kd = 1 ml/g
(b) The case with a constant De

Figure 3.7 Concentration level at CCL base for partially saturated cases with decreasing De .
(Cc = 0.8, k p = 10−9 m/s). Notation: FD: finite deformation model and ND: no deformation
model.
1D finite strain coupled model for consolidation and solute transport 65

3.5.3 EFFECT OF COMPRESSIBILITY OF PORE-WATER (CPW)

As shown in Figure 3.8, the effect of compressibility of pore-water (CPW) is
related to the soil consolidation coefficient. The influence of CPW on the relative
concentration at the CCL becomes more significant for the cases with smaller
consolidation coefficients. When the soil is relatively soft (Cc = 0.8 and k p =
2 × 10−10 m/s), CPW causes twofold longer transit times for the nominal 10%
breakthrough. However, at the early consolidation stage, the retarding effect of CPW
is more pronounced for ‘stiffer’ soils and then the trend reverses (Figure 3.8) after
consolidation completes. These graphs are not shown as the numerical values are
too small to present in the same figure. This can be explained by the slowing fluid
flow and longer consolidation time due to CPW. Since the separation of curves at
a relatively higher concentration level, i.e., absolute concentration difference, is of
interest, it follows that the influence of CPW is more significant in softer soil.
Relative concentration at CCL base, Cf/Cf0

−6
10

−8
10

−10
10
1.5 2 2.5 3
Time, τ (y)

Figure 3.8 Effect of CPW on concentration level at CCL base for partially saturated cases
(Sr = 0.8) with varying De and without sorption (Kd = 0). Solid lines: Cc = 0.8, k p = 2×10−10
m/s; Dashdot lines: Cc = 0.8, k p = 10−9 m/s; Dotted lines: Cc = 0.2, k p = 10−9 m/s. Cross
symbol: with CPW; circle symbol: without CPW (β = 0).

To investigate further the influence of CPW, three models examining the three
terms containing β are considered here.
eSr β ∂ e
• Model A: eliminate (1+e0 )αv ∂t from (3.11);
Sr β e ∂ Q
• Model B: eliminate − 1+e 0 ∂t
from (3.11);
• Model C: eliminate the term involving β from (3.19).

As shown in Figure 3.9, each of the missing terms leads to a large deviation from
the full model, so all terms involving β should be retained for the cases considered.
66 Poro-Elastic Theory with Applications to Transport in Porous Media

0
10

Relative concentration at CCL base, Cf/Cf0

−1
10

−2
10

−3
10
Complete model
Model A
Model B
Model C
−4
10
0 10 20 30 40
Time, τ (y)

Figure 3.9 Significance of each term involving β on concentration level at CCL base for
partially saturated cases (Sr = 0.8, Cc = 0.8, k p = 2 × 10−10 m/s) with varying De and without
sorption (Kd = 0).

3.5.4 EFFECT OF DISPERSION

Lewis et al. (2009) neglected mechanical dispersion on the assumption that the pore
fluid velocity in fine-grain soil is less than 10−6 m2 /s. However, as shown in Figure
3.10, its influence cannot be neglected when the clay is relatively soft, even when
the maximum fluid average linear velocity is approximately 4.5 × 10−9 m/s for the
case Cc = 0.8 and k p = 2 × 10−10 m/s. Its influence becomes more significant as
the hydraulic conductivity increases with the same soil compressibility, Cc . This is
because decreasing De increases the Péclet number (ratio of the rate of advection to
the rate of diffusion). Therefore, a rough estimate using pore fluid velocity alone as
proposed by Lewis et al. (2009) is not always definitive.
Figure 3.11 illustrates the individual influence of decreasing De , dispersion and
CPW. The effect of reducing De causes slower transport, while dispersion a faster
transit. Although the influence of CPW is not as significant as decreasing De and
dispersion, it is not negligible, as shown in Figure 3.11.

3.5.5 EFFECT OF FINITE DEFORMATION

For the soil without sorption (see Figures 3.1(b), 3.7(a), 3.7(b), 3.11(a)), the ND
model always leads to a longer transit time than the finite deformation model. In
the presence of sorption (as shown in Figure 3.7(b)), the difference between the
ND model and the finite deformation model is negligibly small. However, when the
decrease of the effective diffusion coefficient due to deformation is also considered
(Figure 3.7(a) and 3.11(b)), the results of the two models differ.
Compared with the finite deformation model, the small deformation model can
overestimate the contaminant transit time in a liner undergoing large consolidation
1D finite strain coupled model for consolidation and solute transport 67

−1
10

Relative concentration at CCL base, Cf/Cf0

−2
10

−3
10

−4
10
0 5 10 15 20 25 30 35 40
Time, τ (y)

Figure 3.10 Effect of dispersion on concentration level at CCL base for partially saturated
cases (Sr = 0.8) with varying De and without sorption (Kd = 0). Solid lines: Cc = 0.8, k p =
2 × 10−10 m/s; Dashdot lines: Cc = 0.8, k p = 10−9 m/s; Dotted lines: Cc = 0.2, k p = 10−9
m/s. Cross symbol: αL = 0.1 m; circle symbol: αL = 0 (no dispersion).

0
Relative concentration at CCL base, Cf/Cf0

Relative concentration at CCL base, Cf/Cf0

10

−20
10
−1
10

−22
−2 10
10

FD −24
−3
10 FD, CD 10
FD, NLGD
FD, NCPW
−4
ND −26
10 10
0 10 20 30 40 0 0.5 1 1.5 2
Time, τ (y) Time, τ (y)

(a) Kd = 0 (b)Kd = 1 ml/g

Figure 3.11 Comparison of the concentration level at CCL base for various variable
associative in partially saturation soils (Sr = 0.8, Cc = 0.8, k p = 10−9 m/s). Notation: FD: finite
deformation model; CD: constant De ; NLGD: excluding the dispersion; NCPW: excluding the
CPW; ND: no deformation model.

(Figure 3.1(b)). This demonstrates that the significance of geometric nonlinearity

is noticeable for relatively soft soil. This finding is consistent with that of Peters
and Smith (2002) and Lewis et al. (2009). Regarding the consolidation, the small
deformation model can predict settlement that is non-physical for soft soil (i.e.,
larger than the total soil thickness). Therefore, for a relatively compressible soil,
where the consolidation effect is more significant, a finite deformation consolidation
is necessary when being coupled with the solute transport.
68 Poro-Elastic Theory with Applications to Transport in Porous Media

3.6 SUMMARY
In this chapter, a finite deformation model for coupling consolidation and solute
transport processes in partially saturated soil has been presented. It was applied
to predict the VOC breakthrough in a landfill clay liner. CPW, dispersion, the
nonlinear variation of soil compaction, hydraulic conductivity and effective diffusion
are included in the model. Based on the numerical simulation results, the following
conclusions can be drawn.

(1) Consolidation-induced advection has a lasting effect on solute transport

during and after the deformation for relatively compressible soil regardless
of the sorption level, though the sorption can dramatically slow the solute
transport process rate.
(2) After an initial acceleration effect on transport, the finite-deformation
coupled model with decreasing effective diffusion and sorption produces
a lower concentration at the CCL base than the pure diffusion model.
(3) A lower degree of saturation leads to a slower pore fluid flow and solute
transport(since larger pores drain preferentially with decreasing saturation).
The CPW associated with unsaturated conditions cannot be ignored when
the consolidation is required to be coupled with solute transport. In the
model, CPW terms exist in both the consolidation and transport equations,
none of which can be neglected for simplification. Effective diffusion
decreases during consolidation and consequently the relative importance
of mechanical dispersion becomes profound. For a long-term prediction,
mechanical dispersion could cause significant solute transport. Therefore,
it should be included in modeling efforts.
(4) Generally speaking, reducing soil compressibility and improving sorption
levels of clay are the most effective ways to retard contaminant migration.
At the same level of stiffness and sorption, the lower hydraulic conductivity
and lower degree of saturation can lengthen the time for contaminants to
break through the protective liner.