AP08_4.vp


1 Introduction
Old deep mines had to be closed in the neighborhood of

Prague, in the Czech Republic. The problem was that the up-
per part of the cross-section is filled with fresh water and the
lower part with salt water (see Fig. 1). Between these two kinds
of water there is a hygrogeological insulator. During mining
operations, shafts 1000 m deep provided access for miners
were drained, and the salt water was constantly pumped out
of the bottom of the mine. After the mines were closed down,
the shafts also had to be closed to avoid penetration of the salt
water into the upper layer with the fresh water. Therefore,
new stoppers were positioned in the shafts to keep the two
kinds of water separated.

Some parts of the problem to be dealt with are published
recently. Physical modeling of the powder consolidation
processes is published in [1]. This was useful not only for cali-
brating the analytical models but also for identifying the
fundamental mechanisms which form the bases of these mod-
els. A model of the mechanical behavior of the saturated-
-unsaturated porous media is presented within the framework
of the thermodynamics of irreversible processes in [2]. This
model serves as a starting point for formulating the solution
presented in our paper. Some of the formulas for the present
paper are found in [3], where there is an analysis of the math-
ematical model for one-dimensional two-phase flow in porous
media with the aim of finding an exact solution. A generalized
plane strain finite element formulation for the analysis of
poroelastic materials is presented in [4], while the present pa-
per explicitly solves the generalized plain strain in each time
step. Adoptive FE solution strategies allowing for control of
the spatial discretization error of coupled FE-analyses for

porous soils with the pore space filled with water and/or
(compressed) air are given in [5]. In this paper there is a com-
parison of different input data for particular materials. They
served as the basis for the analysis conducted in our paper. In
order to efficiently simulate large-scale contaminant trans-
port flows on serial computing platforms, the authors of [6]
present an explicit finite difference formulation and solution.
Paper [7] addresses partly saturated flow and solute transport
in porous fractured media, and the relevant physical proper-
ties are discussed. Under the assumption of linear soil mass,
nonlinear mechanical behavior or fractions interaction be-
tween groundwater flow and solid deformation is considered.
For a numerical application a sandstone block with and with-
out a fracture is analyzed. The authors of [8] claimed that the
soil-water interaction can significantly influence the forma-
tion and propagation of the bands of localized damage in a
porous medium. In a practically-oriented calculation exam-
ple in [9] shows the need to take into account rock damage,
material softening and dilatation as well as the interaction
between mechanical and hydraulic processes. Paper [10] gives
a brief retrospective review of various analytical methods for
determining the influence on the consolidation process of
vertical drain installations. The results of full-scale tests on
vertical drain behavior are compared with theoretical consoli-
dation rates. Paper [11] describes the change in air water pore
pressures that take place in a freshly compacted soil as the two
phases of the pore fluid stress and solution equilibrium. The
phenomena are demonstrated by a series of measurements
on freshly compacted soils, and are explained by observations
of a simple physical model of the soil. It is concluded that the
processes of equilibration must take place whenever the total
stress on a compacted soil is varied. Paper [12] compares two

34

Acta Polytechnica Vol. 48  No. 4/2008

Modeling Water Pollution of Soil
V. Doležel, P. Procházka, V. Křístek

The government of the Czech Republic decided that in the location to the west of Prague, capital city of the Czech Republic, some deep mines
should be closed because of their low efficiency of coal mined i.e. small amounts and low quality of the coal extracted in the final stage of
mining. The locations near Prague influenced the decision to do maintenance on the abandoned mines, as the thread of soil pollution was
unacceptably high in the neighborhood of the capital city. Before the mines were closed it was necessary to separate existed extensive
horizontal location of salt water below a clay layer in order not to deteriorate the upper fresh water. The salt water could not be allowed to
pollute the upper layer with the fresh water, as many wells in villages in the neighborhood of the former mines would be contaminated. Two
horizontal clay layers (an insulator and a semi-insulator) separated the two horizons containing salt water and fresh water. Before starting
deep mining, vertical shafts had to be constructed with concrete linings to enable the miners to access the depths. The salt water was draining
away throughout the existence of the mine. The drainage was designed very carefully to avoid possible infiltration of salt water into the upper
horizon. Before the mines were abandoned it was necessary to prevent contact between the two kinds of waters in the shafts. Several options
were put forward, the most efficient of which appeared to be one that proposed filling the shafts with spoil soil and creating a joint seal made
of disparate material at the interface between the salt water and fresh water to create a reliable stopper. The material for the spoil soil was
delivered from deposits located not far from the shafts. This material consisted of a variety of grains of sand, big boulders of slate, slaty clay,
sandstone, etc.. Chemical admixtures were considered to improve the flocculation of the filling material. The stopper was positioned at a
depth of 220–300 m below the terrain. As an alternative, thinner stoppers were considered, but this option was discarded.
The aim of this paper is to describe the design of the stoppers applied to separate the two types of water along the contact horizon using Desai’s
DSC theory (Distinct State Concept), and generalized plane strain in the multiphase problem of water flow in a porous medium. In addition,
a comparison of some results from scale experimental models with numerical solutions was carried out. The intrinsic material properties of
stoppers for numerical computations were obtained from physical and chemical laboratory tests. The models were evaluated for the complete
underground work, particularly in its final stage of construction.

Keywords: scale modeling, generalized plane strain, leakage problem.



models of flow in porous media is discussed. The first is the
well known Richard equation, which is based on the assump-
tion that the air in the unsaturated zone has infinite mobility.
This means that it models a single phase. In the second model
a general full two-phase approach is applied, where the air is
considered as a separate phase.

It is worth noting that the coupled (numerical and ex-
perimental) approach can be used for other porous material
insulators and can describe the drying process of such materi-
als (bricks, concrete).

Paper [13] is deals with interaction of subsoil and struc-
tures and includes the influence of pore pressure due to air
and water, and the solid phase. Some results from this paper
are partly used here.

2 Formulation of the problem
In this section we will put forward basic relations leading

to the target formulation of the coupled modeling (experi-
mental and numerical). First, we formulate the problem of the
long-term deformation of the water stopper, which admits ad-
missible seepage, and is denoted as � with the boundary � .
The domain � is equipped with the system of coordinates
0x1x2x3. The problem deals with the distribution of two types
of functions in the domain describing the soil of the stoppers.
The first type of functions are the water pressure p p tw w� ( , )x
and the air pressure p p ta a� ( , )x , which are scalars depend-
ing on the position � �x � x x x1 2 3, , at time t, and the second
type of functions are

� �u x x x x( , ) ( , ), ( , ), ( , )t u t u t u t� 1 2 3 ,

the vector of displacements of the skeleton,

� �u x x x xw w w w( , ) ( , ), ( , ), ( , )t u t u t u t� 1 2 3

is the vector of movements of water in pores and

� �u x x x xa a a a( , ) ( , ), ( , ), ( , )t u t u t u t� 1 2 3

is the vector of movements of the air in the pores. For the two
types of unknowns, two types of equations are available: the
equation of continuity for the water flow in pores and equa-
tions of equilibrium, both being valid at each point x of the
soil. The functions p and u are mutually adjoined, so that the
governing equations create a simultaneous system.

After transforming the coordinate system 0x1x2x3 to the
axisymmetric cylindrical coordinate system 0rz , where r is the
radial coordinate and z is the axial coordinate, the unknown
function becomes:

� �x � r z, , � �u x x x( , ) ( , ), ( , )t u t u tr z� ,

� �u x x xw w w( , ) ( , ), ( , )t u t u tr z� ,

� �u x x xa a a( , ) ( , ), ( , )t u t u tr z� .

The stopper can be considered as a three-phase compos-
ite; this means that the soil contains a solid phase (the skele-
ton), a gaseous phase (the air in the pores), and a liquid phase

(water). The formulation of the problem will first be derived
for the Gaussian orthogonal coordinate system.

The present formulation is restricted to small displace-
ment and small strains. Hence the relationship between the
displacements u and the strains � in the soil skeleton is given
as:

�
	

	

	

	
ij

i

j

j

i

u
x

u

x
� �

�

�

�
�

�

�

�
�

1
2

. (1)

Then for the stresses in the skeleton we can write:


 
 �ij ij ijp� �eff , (2)

where 
ij are the components of the total stress tensor and

ij

eff are the components or the effective stress tensor, both for
i j, , ,�1 2 3in 3D, and�ij is the unit tensor. Relation (1) can for-
mally be written as:

� �� �eff pm, (3)

where the tensors are denoted in vector terms, as usual:

� �� � 
 
 
 � � �11 22 33 12 13 23, , , , , ,

� ��
eff eff eff eff eff eff eff� 
 
 
 � � �11 22 33 12 13 23, , , , , ,

� �m � 1 1 1 0 0 0, , , , , .

Note that the last notation of the components of the
stresses is maintained for tensor notation and is related to the
most common denotation as:

 
x � 11, 
 
y � 22 , 
 
z � 33 , � 
xy � 12 , � 
xz � 13 , � 
yz � 23 .

Now, let us briefly describe the relations being valid for the
deformation of rock mass with the influence of water in pores.
The effective stresses follow the generalized Hooke’s law:

� � � �
eff pl

p� � �D( ) , (4)

where D is the material stiffness matrix of the porous skele-
ton, and �

pl is the plastic or visco-plastic strain tensor. More-
over

�p � �
�

�
��

�

�
��m

p
Kz3

(5)

is the strain of the basic material due to the pore pressure,
and Kz is the bulk modulus of elasticity of the solid material
(matrix).

The total stresses according to (3) to (5) appear to be:

� � � �� � � � �eff plp a pm m( ) , (6)

where

�m
Dm D m m� � � �

�

�
��

�

�
�� 


3
1

3K Kz z
.

The movement of a liquid in the pores is expressed in
term of a time change of the pores
. By superposition of these
influences we obtain in a representative volume element:

� Changes in volume of the skeleton due to �, assuming a
negligible volume compressibility,

� Change in volume (1� n) of the basic material due to p,

35

Acta Polytechnica Vol. 48  No. 4/2008



� Increments of effective stresses cause an increase in the
volume of the skeleton, while the volume of the pores
decreases.
For compression of the liquid due to pore pressure we get

the following relation after easy algebra:


 �� 	 � �
�

� �
�

�
 

!

"
#mT

w�

n
K

n
K

a
K

p
z z

1
1

1
( ) , (7)

where K w is the bulk modulus of elasticity of the liquid and n
is the porosity.

It remains to use the equation of continuity and Darcy’s
law. The equation of continuity expresses the fact that the
amount of liquid is pressed out of the pores in the deformed
soil in a unit of time. This has to be equal to the increase in the
volume of pores in the same time unit (due to the law of mass
conservation):

ndivv �� �
, (8)

where v is the vector of the relative velocity of the displace-
ment of the liquid with respect to the displacement of the
skeleton.

36

Acta Polytechnica Vol. 48  No. 4/2008

Fig. 1: Cross-section of the shafts and insulators
– original state after abandoning the mine
– geological conditions in the overburden of the mine



According to Darcy’s law the relative velocity of the flowing
liquid is proportional to the gradient of the pore pressure:

n
g
k

pf f f f� �v g�� �grad , (9)

where k is the coefficient of filtration (permeability),

� �g � 0 0, , g is the vector of gravitation acceleration, providing
the third coordinate is vertical, � f is the density of the phase,
f w a� , and n is the porosity. The last three relations lead us
to the following equation:

� � � �
c

K
p ps � � �� �mT

� 0, (10)

where

� �� �
�

� �
n

K
n

K Kz z
w

1
1

1
( ) ,

and

c
K k

g

s

k
�

�

is the coefficient of consolidation, and K s is the bulk modulus
of elasticity of the skeleton. The overdot means differentia-
tion by the time, and g is the gravitation acceleration.

In the sense of the assumption stated above that there is a
consolidated stage let us set � �0. The final governing equa-
tion for the liquid pressure in pores is written as:

� � �
c

K
p ps � � � 0. (11)

The latter equation (11) has to be rearranged in such a
way that different properties are considered in different di-
rections. Then the steady state equation (11) turns into the
following equation:

� �

�
�

�

�
� �

�

�
��

�

�
�� � �

�
�

�

�
� �

	

	

	

	

	

	

	

	

	

	

	

	

	

	x
k

x y
k

y z
k

zx y z t
p

$
%
&

'
(
)

�0, (12)

when assuming ki, i x y z� , , , to be dependent on the position
of the point x in the domain � . If ki is different from kj for at
least one case i j* then the problem is anisotropic.

If we take into consideration saturation of the phases, the
effective stresses are expressed as:

� �
eff w w a a� � � �( ) ( )1 n n S p S p m, (13)

where Sw is the degree of saturation of the water phase and Sa

is the degree of saturation of the air. It obviously holds:
S Sw a� �1. Considering this, equation (4) can be expressed
as:

	


	
�

ij

jx

eff

� �g 0, (14)

where

� �� � � �( )1 n nS p nS ps w w a a

is the average density of the three-phase medium and g is the
gravity acceleration vector. The superscripts refer to the solid
phase s, the water phase w and the air phase a.

Darcy´s law for the artificial velocity �v f of a fluid f w a� ,
relative to the soil skeleton can be derived from the linear
momentum balance equation for the respective fluid phase
as, see [5]:

� ��v gf
f

f
f fk

p g
p p� �( )grad T (15)

with k f denoting the coefficient of permeability. If two fluid
phases (either water or air) are present in the pores, then k f

depends on the degree of saturation S f .

Considering a domain without any disturbances, the mass
balance equation for a fluid phase f yields, [5]:

� � � �f f f f f f f fS nS n S� � � ( )�vol div� � � �
�v (16)

or alternatively:

	

	
� �

t
nS f f f f( ) ( )� �div �v . (16a)

In (16) the term on the right hand side represents the in-
flow of fluid mass into a given volume element which has to be
balanced by the terms on the left hand side of this equation;
i.e. the inflow can be stored in the volume element by an in-
crease in the volumetric strain of the soil skeleton ��vol by an
increase in the density �� f and by the degree of saturation �S f

of the respective fluid phase f. The constitutive equation for
the compressible barotropic fluid phases is given as:

� ��

�

f

f

f

f
p

K
� � (17)

37

Acta Polytechnica Vol. 48  No. 4/2008

Fig. 2: Alternative setting of more than one stopper



with K f the bulk modulus of the respective fluid phase f. The
minus sign on the right hand side of (17) is a consequence of
the definition of tensile stress as a positive quantity. Assuming
the air phase at constant temperature it follows from (17) us-
ing Boyle-Marriote’s law pV � const. (V is the volume of � )
and the principle of conservation of mass, that the bulk
modulus K a is equal to the magnitude of the air pressure pa.

Summing up the above theory, the first system of equa-
tions can be formulated as:

� �	

	

 �

x
S p S p b

j
ij ij i� � � � �( )w w a a 0, (18)

where

� �� � � � �
 
 �ij ij ijn S p S p( ) ( )1 w w a a ,

�ij is again the unit tensor of the second order and

� �b n nS nS gi � � � �( )1 � � �s w w a a , see (14).

The boundary conditions of � are defined by either the
given displacements u or the given tractions t .

Substituting (15) to (16) yields completion of the system of
equations:

1
�

	

	

	

	

	

	

	

	

	

	

	

	f x
f

y
f

z
f

g x
k

x y
k

y z
k

z
� �

�
�

�

�
� �

�

�
��

�

�
�� � �

�
�

�

�
�

�

�
 

!

"
#

� � � �

p

n
S
K

p S f w a
f

f
f f� � , ,�vol 0

(19)

or alternatively

1
�

	

	

	

	

	

	

	

	

	

	

	

	f x
f

y
f

z
f

g x
k

x y
k

y z
k

z
� �

�
�

�

�
� �

�

�
��

�

�
�� � �

�
�

�

�
�

�

�
 

!

"
#

� � � �

p

n
S
K

p S n f w a
f

f
f f� � , , .0

(19a)

The boundary conditions are either the given hydrostatic
fluid stress or the averaged flux of the respective fluid phase
through the boundary of the domain � .

Comparing (19) with (12) yields that if the volumetric
change of strains is zero and the increment of saturation is
also zero, equation (12) is equal to (19) for the condition given
by:

� f f

f
g nS

K
�1. (20)

3 Weak formulation
Now we very briefly derive the weak formulation of the

problem (12) extended by the nonhomogeneous right hand
side f (which involved influence of saturation, porosity, etc.)
from which the variational formulation follows. The bound-
ary conditions wil1 be a consequence of this formulation.
Multiplying the last equation by a smooth enough function �,
which is not identically equal to zero in the domain �, and
after integrating the derived equation over the domain we ar-
rive at:

�
�

�
�

�

�
� �

�

�
��

�

�
�� �

�

�
�

�

�
� �

	

	

	

	

	

	

	

	

	

	

	

	

	

x
k

p
x y

k
p
y z

k
p
zz

p
t

f

	
�

�

$
%
&

'
(
)

�

+

+

d

d

�

�

�

�

,

(21)

Using Green’s theorem the latter equation turns to:

k
p
x x

k
p
y y

k
p
z z

p
t

p
z

	 	�

	 	

	 	�

	 	

	 	�

	 	

	

	
�

	

	
� � �

$
%
&

'
(
)

�+ d�

�
n

f

d�

�

�

�

+

+� �d ,

(22)

where
	

	

	

	

	

	

	

	n
n k

x
n k

y
n k

zx y z z� � �

and � �n n n nz y z� , , is the outward unit normal to the bound-
ary � of � . From (13) the weak formulation follows and also
the finite element method starts with this equation. The
boundary conditions can also be derived as:
	

	

p
n

pn pis prescribed on and is prescribed on� � ,

,� � � � �n p n p, � � -0
(23)

(since the permeability coefficients are also defined by vector
transformation) after some algebra we eventually get:

r
r

k
p
r r

k
p
r z

k
p
z

p
z�

�

�
�

�

�
� �

�

�
�

�

�
� �

�

�
�

�

�
� �

	

	

	

	

	

	

	

	

	

	

	

	

1
t

rf

$
%
&

'
(
)

�

� �

�

�

+

+

�

�

d

d

�

�

�

�

,

(24)

which is valid for arbitrarily � from an admissible set of solu-
tions. Such an admissible set can be created from continuous
functions, which fulfill the boundary conditions (23). The do-
main �� is now a cube (0, R) × (0, 2� ) × (0, H), where H is the
height of the water stopper, and R is the radius.

As the problem is considered axisymmetrical, neither the
given quantities nor the searched quantities depend on the
circumferential coordinate. Consequently, the cube turns into
a rectangle. Dividing the rectangle into finite elements of rect-
angular shape, a linear distribution of the function p and con-
stant distribution of

	

	

p
n

is obtained under the condition that
p is bilinear oon each rectangle. Introducing this approxima-
tion, the folloving equations are obtained:

A p p f� �
	

	t
(25)

where p is the vector of nodal pressures and f involves the
influence of changes of internal characteristics during the
process of settlement inside the stopper.

The last equation is algebraic simultaneous equation ex-
cept for the time dependent vector. The latter term may be
treated in the sense of a finite difference scheme. Consider

38

Acta Polytechnica Vol. 48  No. 4/2008



one typical time interval �t t ti i i� � �1 and suppose that the
solution pi�1 is known at time ti�1. Inside this interval we apply
a linear approximation for the vector p p� ( )t :

p p p( ) ( ) ( ) ( )t t ti i� � � �� �1 1 (26)

where � � � �( )t t ti 1 � . In the same way we approximate the
vector f. In the sense of the above introduced approximations
the time derivative takes the form:

	

	

p p p
t t

i i�
� �1
�

. (27)

Substituting (23) to (21) yields:

� �( ) ( )( )A p Bq
p p

A Bqi i
i i

i it
� � �

�
� � ��

� �
1

1 11
�

p . (28)

4 Scale modeling
The mathematical model serves for solving water pres-

sures in pores and the overall (total) stresses in the insulator.
Note that the components of the stress tensor are virtually
impossible to determine either in experimental or “in situ”
measurement. On the other hand, no numerical model can
be applied without appropriate material and internal pa-
rameters. For this reason experimental scale models were
prepared and the process of constructing the insulator was
simulated in the stands, which are partly glassed basins boxes.
The sides of the stands are about 1.5×1.5×1.5 m2. Natural
materials fill the stands and the water is pushed through the
pores of the material.

Classification of the insulator material originates from the
description of the samples from the site and from laboratory
tests. Due to unexpected inflow of ground water and in order
to save energy, continuation of drainage was excluded and the
project assumed no penetration of saline water into the upper
level. The problem consists of the time dependent behavior
of the insulator, particularly its sedimentation, compression,
seepage, change of coefficients of filtration (k k kt z� �� 08. ; all
these coefficients are dependent on z, which is axial – verti-
cal – coordinate), and the solution is strongly nonlinear. All
these factors essentially influence the insulation function of
the water stopper. Under the assumption of decreasing inflow
of water into the shafts in dependence on the increasing water
level, and the velocity of construction of the stopper (20 to
40 m a day) it is realistic to imagine such a situation that the
salt water reaches the upper level of the stopper at some stage
of construction of the stopper. Neither is this case nor if
the stopper is finished, the two types of water should not
come into contact, i.e., the upper level of saline water and the
lower level of the fresh water must not meet at any stage of
construction.

From extensive experimental measurements on scale
models aimed at simulating the stoppers in the shafts a very
important decision follows, which is often mentioned in pub-
lications and discusses the validity of classical form of Darcy’s
law, see Fig. 3.

It was experimentally proved on physical scale models
that during leakage through a fine-grained material with low
permeability there is a hygraulic gradient limit value at which
the flow stops. This is caused by great intermolecular forces
that adhere the water to the walls of pores. This means that

there is a flow after a limit value of the hygraulic gradient is
exceeded. For the materials of the stoppers this value ranges
between 0.8 and 0.95 in our case. The second restriction in
the applicability of Darcy’s law follows from leakage thought
very course-grained material, which causes turbulent flow.

Reynold’s number is used to distinguish these two states.
In very fine-grained, clayey materials (which were supposed
to be for use in the stoppers) laminar flow occurs if Reynold’s
number ranges between 2 and 8. The viscous forces are preva-
lently inertial in this case.

To ensure an appropriate material for the stoppers, a spe-
cial technology had to be developed. For such a technology,
an electrochemical couple had to be provided while successive
backfill sediments were created. It appeared that the engi-
neering has to obey the law that the so-called salt flocculation
is prevalent. That is, between each two particles of the clayey
backfill the attractive forces are greater than the repulsive
forces (disperse structure). The process of flocculation de-
pends on double electric layers between distinctive clay parti-
cles. From the laboratory tests on the clay materials it was
derived that the moisture required for the stoppers has to be
lower than the flow index WL. The flow index ranges from
35.9 % to 81.4 % for the given materials. The flow index can
be determined from the relationship of the moisture and the
coefficient of permeability k, which is shown in Fig. 4. The
graph in Fig. 4 holds valid for materials with a flow index
approaching 80 %. A subsequent drop in moisture and con-
sequent increase in the impermeability of the stopper can be
attained by flocculation admixtures, which speed up the pro-
cess of sedimentation and contribute to the creation of a
stiffer linkage between distinctive particles. These conclusions
were confirmed by the results of laboratory tests performed
on site, which aimed to state the sedimentation and floccula-
tion regimes for various types of rock and soil used for the
stopper.

This appears to be very important for the identifying the
material relation between the effective porosity and the aver-
age size of the grains, Fig. 4. This relationship determines the
design of the construction technology of the stoppers.

39

Acta Polytechnica Vol. 48  No. 4/2008

Fig. 3: Dependency of low velocity on hydraulic gradient



Fig. 5 depicts the assumed time dependency on the poros-
ity. This picture shows that time limits are reached at the level
of porosity required for restricting mixture of salt water and
still water. Note that the porosity decreases due to floccula-
tion, and can be speeded up if sufficient money is made
available.

Scale modeling proves that correct functioning of the
stopper is very tough and difficult to achieve. Moreover, scale

modeling proved the material characteristics changed signifi-
cantly in the final stage, e.g., Young’s modulus E dropped by
50 to 60 percent and toughness by 30 to 40 percent.

5 Desai’s model
Desai’s model is introduced because of the different vol-

ume fractions of the pores and skeleton in different phases of
the process of constructing the stopper. It is not possible to
feed an algorithm of computation by reasonable input data to
get unique results; the problem is strongly nonlinear, as stated
above, so that also input information sensitive.

Consider an element of porous (soil) material with a fluid
(water). In the theory of composite materials we say that RVE
(the reference volume element) is considered. The applied or
total compressive force, P, is carried by the force in the solid
skeleton, F s, and is transmitted through the contact area, As.
Force F f is transmitted through the fluid area Af. Then the
force equilibrium gives:

P F F� �f s . (29)
Dividing the last relation over the nominal area, A, of the

RVE, after rearrangement we obtain:
F
A

F
A

A
A

F
A

A
A

a f

f

f s

s

s
� � , (30)

where A
A

f
and A

A

s
are volume fractions of fluid and solid.

This denotation from Desai’s composite theory, [14], substi-
tutes as:

� � �
a f s� � �( )1 D D . (31)

The latter denotation is advantageous from certain points
of view. Since the volume fractions sometimes change basi-

40

Acta Polytechnica Vol. 48  No. 4/2008

Fig. 4: Dependency of effective porosity on average grain size

Fig. 5: Dependency of time on porosity

Fig. 6: Typical distribution of the disturbance function

middle grain sandstone sand siltstone silt sandstone clay siltstone

shaft dump shaft dump shaft dump shaft dump

Rd 46.9 35.2 78.3 15.1 67.6 66.1 15.1 15.1

Rtp 2.47 2.98 2.9

E 7487 10267 2550 18063 6280 9528

E0 6603 1200 10209 14816 7366

v 0.19 0.26 0.36 0.4

Rsh 4.9 6.7 9.9 7.55

� 61.6 28 65.2 31.2 60.2 67.9

Table 1: Material properties of the spoil soil



cally, function D can describe the real situation more naturally
than in the case of classical composites. This function is called
the disturbance function and can be determined from scale
or “in situ” measurements from the condition that the (lin-
early, say) computed results in some selected points and the
measured results should be as close as possible at the same
points. A typical distribution of the disturbance function in
our problem is depicted in Fig. 6.

Assume, for example, that only water pressure on the
stopper is considered on its upper boundary. Then (31) can
be written as:


 
 � 
z z zD h Da f s� � � �( )( ),1 0 , (32)

where h is the distance from the upper level of the stopper
downwards, � is volume weight and
z,0

f is the pressure on the
upper boundary of the insulator.

6 Some results and discussion
We first introduce the very important distribution of coef-

ficients of filtration k (m/s) in the vertical direction depending
on moisture w (%), which was derived from large-scale mod-
els. The relation is depicted in Fig. 7.

The assumptions made in the numerical test models are in
compliance with the experimental data and are used, first, to
improve the model, and, second, as input data for feeding the
examples.

The coefficients of saturation of water and of air are con-
sidered as: Sw �092. , Sa �008. . Note that the salt water must
be reflected as artesian, i.e. it is pushed to the stopper from
below, so that air pressure also appears in the formulas.

The bulk modules are listed as:

K w � .215 106. kN/m3,
K a � .10 102. kN/m3,
Kz � .485 106. kN/m3.

The densities of the constituencies are given by:

�w �10 kN/m3,
�a �00214. kN/m3,
� s �170 kN/m3.

The equations that have to be solved are listed below, with
the coefficients mentioned above:

The first equations read:

� �	

	

 �

x
S p S p b

j
ij

w w
ij i� � � � �( )a a 0

where

� �� � � � �
 
 �ij ij
w w

ijn S p S p( ) ( )1 a a ,

�ij is again the unit tensor of the second order and

b n nS nS gi
w w� � � �[( ) ]1 � � �s a a .

The second system of equations is:

1
�

	

	

	

	

	

	

	

	

	

	

	

	f x
f

y
f

z
f

g x
k

x y
k

y z
k

z
�

�

�
��

�

�
�� �

�

�
��

�

�
�� �

�

�
��

�

�
��

�

�
 

!

"
#

� � � �

p

n
S
K

p
t

S n f w a
f

f
f f� ( ) , , .

	

	
0

The distribution of water and air pressures, and the
stresses from Lame’s equations are computed by finite ele-
ment method. At certain points, strength in the sense of the
Mohr-Coulomb hypotheses (ideally elasto-plastic law) is at-
tained and the disturbance function requires a change in its
values at such points (having an impact on the coefficients of
filtration and others). The resulting water pressure of salt wa-
ter in the vertical and radial directions appears in Table 2.
Note that the radius of the shafts is R � 4.45 m. It is clear that
the highest velocity values are observed at the outer ring of
the stopper and on the lower boundary. Moreover, 16 m from
the lower boundary almost no pressure is obtained in the final
stage of the construction of the stopper, the distance in z
direction is measured in m. Radius r is also measured in
distance in meters from the shaft axis. The distribution of
the water pressure along the lower boundary shows that the
highest pressure is observed at the concrete lining, and the
lowest at the axis.

7 Conclusions
Based on the results of experiments and numerical analy-

sis, the following conclusions can be formulated:

� The optimal value of the coefficient of permeability should
range between 10�9–10�11 m/s. Such materials have been
prepared for use in sufficient quantities.

� From the point of view of the future insulation properties
of the stoppers the backfill should be carried out in a dry
space in the shaft.

� A possible negative effect of the friction along the walls of
the shaft and the backfill material could be increased by a
drop in the angle of integral friction, i.e. by increasing the
moisture, adding bentonite, etc. However, the real situa-
tion on site requires such a technology for constructing the
stoppers such that the backfill material is deposited in a
very high water column (400�600 m).

� Engineering of the stoppers with the required properties is
strictly dependent on the moisture of the material, which
must not exceed a certain limit value in the final stage of
constructing the stoppers.

41

Acta Polytechnica Vol. 48  No. 4/2008

0,00

2,00

4,00

6,00

8,00

10,00

12,00

14,00

16,00

18,00

20,00

0 10 20 30 40 50 60 70 80 90 100 110

W

x
.1

0
e

x
p

(-
1

1
)

Fig. 7: Relation coefficient of filtration k and moisture w



� The moisture value should in every case be less than the
limit of flow and greater than the limit of plasticity. This
depends on the sedimentation velocity and on the creation
of a stiff structural linkage between the particles. This
process can be very efficiently supported by flocculant’s
admixtures.

� A comparison of the pressure value at different cross sec-
tions of the stopper yields that when the thickness of the
stopper is increased above 80 m the permeability does not
increase, as the entire amount of flow disappears to the
sides of the stoppers.

� From the mathematical solution using FEM it follows that
for a stopper thickness of approximately 80 m and for a
hygrostatic pressure below 80 MPa, water flow amounts to
0.711 l/day for unit porosity.

� If more stoppers are added, their insulating effect is negli-
gible.

Appendix

Solution of the problem of
permeability of the stoppers in
cylindrical coordinates
Since the problem should involve rotationally symmetric

anisotropy (the original concrete lining has to appear in the
formulation of the problem), a transformation to cylindrical
coordinates leads to more appropriate expressions. Consider
first the 2D harmonic anisotropic operator and introduce
polar coordinates. The cylindrical coordinates will then be

supplemented by natural addition of the third direction.
Hence, the following operator

�2

2

2

2

2� �
	

	

	

	x y
(33)

has to be changed by introducing a transformation
x r� cos � , y r� sin� . (34)

It obviously holds
	

	

	

	

	

	

	�

	

	

	�x
r
x r x

� � , 	

	

	

	

	

	

	�

	

	

	�y
r
y r y

� � , (35)

or
	

	

	

	

	

	

	

	

	

	

	

	

	

	x

y

r
x x
r
y y

r
$

%
//

&
/
/

'

(
//

)
/
/

�

�

�

 
 
 
 

!

"

#
#
#
#

	

	�

$

%
//

&
/
/

'

(
//

)
/
/

. (36)

Inverse formulation to (21) gives:
	

	

	

	

	

	

	

	

	

	r
x
r x

y
r y

� � , 	

	�

	

	�

	

	

	

	�

	

	
� �

x
x

y
y

or
	

	

	

	�

	

	

	

	

	

	�

	

	�

r
x
r

y
r

x y

$

%
//

&
/
/

'

(
//

)
/
/

�

$

%
//

&
/
/

'

(
//

)
/
/

$

%
//

&
/
/

'

(
//

)
/
/

�

$

%
//

&
/
/

'

(
//

)
/
/

	

	

	

	�

	

	

	

	

r x

y

T (37)

where

T �
�

�

�
 

!

"
#

cos sin
sin cos

� �

� �r r
T =

Since T T� �1 T, (26) can be rewritten as:

42

Acta Polytechnica Vol. 48  No. 4/2008

z/r 0.01 0.6 1.2 1.8 2.4 3.2 3.8 4.45

16 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

15 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001

14 0.000 0.000 0.001 0.002 0.002 0.001 0.001 0.002

13 0.002 0.002 0.002 0.004 0.004 0.004 0.005 0.005

12 0.004 0.007 0.011 0.013 0.015 0.015 0.017 0.019

11 0.014 0.021 0.034 0.046 0.054 0.057 0.057 0.057

10 0.037 0.060 0.107 0.148 0.175 0.189 0.194 0.195

9 0.080 0.121 0.205 0.279 0.329 0.354 0.365 0.369

8 0.150 0.223 0.378 0.518 0.609 0.658 0.673 0.675

7 0.271 0.406 0.690 0.950 1.122 1.208 1.236 1.239

6 0.483 0.725 1.248 1.733 2.056 2.221 2.266 2.265

5 0.835 1.262 2.212 3.127 3.758 4.083 4.159 4.148

4 1.358 2.084 3.759 5.496 6.801 7.518 7.684 7.650

3 2.004 3.131 5.877 9.062 11.908 13.837 14.461 14.438

2 2.542 4.050 7.931 13.020 18.778 24.742 28.426 29.187

1 2.597 4.205 8.440 14.416 22.530 35.209 58.596 74.068

Table 2: Salt water pressure at different points of the stopper



	

	
�

	

	

� 	

	�x r r
� �cos

sin , 	

	
�

	

	

� 	

	�y r r
� �sin

cos (38)

or

	

	

	

	

	

	

	

	�

x

y

r
$

%
//

&
/
/

'

(
//

)
/
/

�

$

%
//

&
/
/

'

(
//

)
/
/

�T 1 .

Using (23) the harmonic operator is expressed as:

� �

$

%
//

&
/
/

'

(
//

)
/
/

$

%
//

&
/
/

'

(
//

)
/
/

�

	

	

	

	

	

	

	

	

	

	

	

	�

x

y

x

y

r

T
$

%
//

&
/
/

'

(
//

)
/
/

$

%
//

&
/
/

'

(
//

)
/
/

� �

T

TT T

	

	

	

	�

	

	

	

	

r

r r r

2

2
1

�
	

	 �

2

2 .

(39)

So far, an isotropic harmonic operator has been consid-
ered. For an anisotropic harmonic operator a similar ap-
proach can be applied, but natural transformation of the coef-
ficients of anisotropic permeability

k k kx r� �cos sin� �� , k k ky r� �sin cos� ��

is not used, as the transformation itself covers the coefficient
of rotational anisotropy. Substituting the above formulas to
the anisotropic harmonic operator in 2D yields

�2 p
x

k x y
p
r r

p

y

r� �
�

�
��

�

�
��

�

�
 

!

"
#

�

	

	
�

	

	

� 	

	�

	

	

( , ) cos
sin

k x y
p
r r

p

r

� �
	

	

� 	

	�

�
	

	

( , ) sin
cos

cos
si

�
�

�
��

�

�
��

�

�
 

!

"
#

� �
n

( , ) cos
sin

si

� 	

	�
�

	

	

� 	

	�r
k x y

p
r r

p
r

�

�
��

�

�
�� �

�

�
��

�

�
��

� n
cos

( , ) sin
cos

�
	

	

� 	

	�
�

	

	

� 	

	�
�r r

k x y
p
r r

p
�

�

�
��

�

�
�� �

�

�
��

�

�
��

� �

�

cos
sin cos

sin cos

2 �
	

	

	

	

� � 	

	�

	

	

� � 	

r
k

r
p

r
k

r
p

r

r r

	

	

	�

� 	

	�

	

	�

�
	

	

� �

r
k p

r
k p

k
p

r r

r r

r

�

� �

sin

cos
sin cos

2

2

2
2

2 2
	

	 	�

� 	

	 �

�
	

	

	

	

�
�

2 2

2

2

2

2

p
r r

p

r
k

r
p

�
�

�
�
�

�

�
�
�

� �

sin

sin
sin cos

sin cos cos

� 	

	

	

	�

� � 	

	�

	

	

� 	

	�

	

	�

�

�

r r
k p

r
k

r
p

r
k py� �

�

2

2

k
p

r r
p

r r

p
y sin

sin cos cos2
2

2

2 2

2

2

22�
	

	

� � 	

	 	�

� 	

	 �
� �

�

�
�
�

�

�
�
�.

(40)

In the isotropic case, k k k k kx y r� � � �� , depending on
the position and the time. k� and kr are the permeability coef-
ficients in the hoop and radial directions, respectively.

In the case of axisymmetry with respect to the axis of the
shaft we obtain in cylindrical coordinates in the sense of
distributions:

�2

2

2p r z t k r z
p r z t

r

k r z
r

p r z t
r

r( , , ) ( , )
( , , ) ( , ) ( , ,

� �
	

	

	

	

	 )

( , )
( , , ) ( , ) ( , , )

.

	

	

	

	

	

	

	

r

k r z
p r z t

z

k r z
z

p r z t
zz

z� �
2

2

(41)

Acknowledgment
This paper was prepared with financial support

from GAČR 103/05/0679, and from GAČR project No.
103/05/2227. The second author is also sponsored by a grant
of the Ministry of Education of the Czech Republic number
MSM6840770001.

References
[1] Piehler, H. R., DeLo, D. P.: Physical Modeling of Pow-

der Consolidation Processes. Progress in Materials Science,
Vol. 42 (1997), p. 263–276, UK.

[2] Schrafler, B. A.: Multiphase Flow in Deforming Porous
Material. International Journal for Numerical Methods in
Engineering, Vol. 60 (2004), p. 27–50.

[3] Fučík, R., Beneš, M.: Numerical Analysis of Multiphase
Filtration Flow in the Problem of Underground Water
Pollution (in Czech: Numerická analýza vícefázového
filtračního proudění v problematice znečištění podzem-
ních vod). Research project, CTU, Faculty of Nuclear and
Physical Engineering, Prague 2005.

[4] Cui, L., Kaliakin, V. N., Abousleiman, Y., Cheng A. H.
-D.: Finite Element Formulation and Application of
Poroelastic Generalized Plane Strain Problems. Inter-
national Journal of Rock Mechanics and Mining Science,
Vol. 34 (1997), No. 6., p. 953–962, UK.

[5] Stelzer, R., Hofstetter, G.: Adaptive Finite Element Anal-
ysis of Multi-Phase Problems in Geotechnics. Geotechnics
and Geomechanics, Vol. 32 (2005), p. 458–481, UK.

[6] Rao, P., Medina, M.: A Multiple Domain Algorithm
for Modeling Two Dimensional Contaminant Trans-
port Flows. Applied Mathematics and Computation, Vol. 174
(2006), Issue 1, p.117–133, USA.

[7] Schmid, G., Rosen, B.: Coupled Ground-Water Pro-
cesses in Porous Fractured Media. Computer Method and
Advances in Geomechanics, (2001), p. 1755–1760, A. A.
Balkema, Rotterdam.

[8] Liu, X., Scarpas, A., Blaauwendraad, J.: Finite Element
Investigation of Localization in Porous Media. Computer
Method and Advances in Geomechanics, (2001), p. 559–564,
A. A. Balkema, Rotterdam.

[9] Hou, Z., Lux, K.,-H.: Analysis of Hydro-Mechanical
Interaction of Rock Salt in Underground Waste Reposi-
tories. Computer Method and Advances in Geomechanics,
(2001), p. 647–652, A. A.Balkema, Rotterdam.

[10] Hansbo, S.: Design of Vertical Drain Installations Based
on Darcian and Non-Darcian Flow: A Comparison of
Real Behavior and Theory. In: GEOTECH-YEAR 2000,
Developments in Geotechnical Engineering, (2000)
p. 179–192, Bangkok, Thailand.

[11] Blight, G. E.: Air-Water Solution Processes in Recently
Compacted Soil. GEOTECH-YEAR 2000, Developments

43

Acta Polytechnica Vol. 48  No. 4/2008



in Geotechnical Engineering, (200), p. 303–310, Bangkok,
Thailand.

[12] Tegnander, C.: Models for Ground Water Flow: A Nu-
merical Comparison Between Richards’ Model and the
Fraction Flow Model. Transport in Porous Media, Vol. 43
(2001), p. 213–224, Netherlands.

[13] Krejčí, T., Nový, T., Sehnoutek, L., Šejnoha, J.: Struc-
ture – Subsoil Interaction in View of Transport Processes
in Porous Media. CTU Report, 1, Vol. 5 (2001).

[14] Desai, Ch.: Mechanics of Materials and Interfaces: The Dis-
turbed State Concept. CRC Press LLC 2001.

Doc. Ing. Vladimír Doležel, CSc.
e-mail: vladimir.dolezel@upce.cz

Geotechnics and Foundation Engineering, Transport
Infrastructure Reconstruction and Maintenance

Transport Faculty of Jan Perner
University of Pardubice
Studentská 95
532 10 Pardubice, Czech Republic

Prof. Ing. RNDr. Petr Procházka,DrSc.
e-mail: petrp@fsv.cvut.cz

Department of Mechanics

Prof. Ing. Vladimír Křístek, DrSc.
e-mail: kristek@fsv.cvut.cz

Department of Concrete and Masonry Structures

Czech Technical University in Prague
Faculty of Civil Engineering
Thákurova 7
166 29 Praha 6, Czech Republic

44

Acta Polytechnica Vol. 48  No. 4/2008