Modeling Water Pollution of Soil

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.


Introduction
Old deep mines had to be closed in the neighborhood of Prague, in the Czech Republic.The problem was that the upper 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 calibrating the analytical models but also for identifying the fundamental mechanisms which form the bases of these models.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 mathematical 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 paper 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 comparison 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 transport 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 properties are discussed.Under the assumption of linear soil mass, nonlinear mechanical behavior or fractions interaction between groundwater flow and solid deformation is considered.For a numerical application a sandstone block with and without a fracture is analyzed.The authors of [8] claimed that the soil-water interaction can significantly influence the formation and propagation of the bands of localized damage in a porous medium.In a practically-oriented calculation example 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 consolidation 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

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 assumption 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 experimental) approach can be used for other porous material insulators and can describe the drying process of such materials (bricks, concrete).
Paper [13] is deals with interaction of subsoil and structures 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.

Formulation of the problem
In this section we will put forward basic relations leading to the target formulation of the coupled modeling (experimental and numerical).First, we formulate the problem of the long-term deformation of the water stopper, which admits admissible seepage, and is denoted as W with the boundary G.The domain W is equipped with the system of coordinates 0x 1 x 2 x 3 .The problem deals with the distribution of two types of functions in the domain describing the soil of the stoppers.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 equations 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 0x 1 x 2 x 3 to the axisymmetric cylindrical coordinate system 0rz , where r is the radial coordinate and z is the axial coordinate, the unknown function becomes: = .
The stopper can be considered as a three-phase compos-ite; this means that the soil contains a solid phase (the skeleton), 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 displacement and small strains.Hence the relationship between the displacements u and the strains e in the soil skeleton is given as: Then for the stresses in the skeleton we can write: where s ij are the components of the total stress tensor and s ij eff are the components or the effective stress tensor, both for i j , , , =1 2 3in 3D, and d ij is the unit tensor.Relation (1) can formally be written as: where the tensors are denoted in vector terms, as usual: { } s = s s s t t t 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: 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: where D is the material stiffness matrix of the porous skeleton, and e pl is the plastic or visco-plastic strain tensor.Moreover is the strain of the basic material due to the pore pressure, and K z is the bulk modulus of elasticity of the solid material (matrix).
The total stresses according to (3) to (5) appear to be: where .
The movement of a liquid in the pores is expressed in term of a time change of the poresq.By superposition of these influences we obtain in a representative volume element: • Changes in volume of the skeleton due to e, assuming a negligible volume compressibility, • Change in volume (1n) of the basic material due to p, • 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: 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): where v is the vector of the relative velocity of the displacement of the liquid with respect to the displacement of the skeleton.According to Darcy's law the relative velocity of the flowing liquid is proportional to the gradient of the pore pressure: where k is the coefficient of filtration (permeability), { } g = 0 0 , , g is the vector of gravitation acceleration, providing the third coordinate is vertical, r 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: where and is the coefficient of consolidation, and K s is the bulk modulus of elasticity of the skeleton.The overdot means differentiation 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 a » 0. The final governing equation for the liquid pressure in pores is written as: The latter equation (11) has to be rearranged in such a way that different properties are considered in different directions.Then the steady state equation (11) turns into the following equation: when assuming k i , i x y z = , , , to be dependent on the position of the point x in the domain W. If k i is different from k j 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: where S w is the degree of saturation of the water phase and S a is the degree of saturation of the air.It obviously holds: . Considering this, equation ( 4) can be expressed where 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 relative to the soil skeleton can be derived from the linear momentum balance equation for the respective fluid phase as, see [5]: [ ] Considering a domain without any disturbances, the mass balance equation for a fluid phase f yields, [5]:  Summing up the above theory, the first system of equations can be formulated as: where , d ij is again the unit tensor of the second order and , see (14).
The boundary conditions of G are defined by either the given displacements u or the given tractions t .
Substituting (15) to ( 16) yields completion of the system of equations: 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 G.
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:

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 boundary conditions wil1 be a consequence of this formulation.Multiplying the last equation by a smooth enough function j, which is not identically equal to zero in the domain W, and after integrating the derived equation over the domain we arrive at: Using Green's theorem the latter equation turns to: where is the outward unit normal to the boundary G of W . From ( 13) the weak formulation follows and also the finite element method starts with this equation.The boundary conditions can also be derived as: (since the permeability coefficients are also defined by vector transformation) after some algebra we eventually get: which is valid for arbitrarily j from an admissible set of solutions.Such an admissible set can be created from continuous functions, which fulfill the boundary conditions (23).The domain ¢ W is now a cube (0, R) × (0, 2p ) × (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 rectangular shape, a linear distribution of the function p and constant distribution of ¶ ¶ p n is obtained under the condition that p is bilinear oon each rectangle.Introducing this approximation, the folloving equations are obtained: 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 except for the time dependent vector.The latter term may be treated in the sense of a finite difference scheme.Consider one typical time interval Dt t t i i i = --1 and suppose that the solution p i-1 is known at time t i-1 .Inside this interval we apply a linear approximation for the vector p p = ( ) t : where l = -- ( ) t t t i 1 D .In the same way we approximate the vector f.In the sense of the above introduced approximations the time derivative takes the form: Substituting ( 23) to (21) yields:

Scale modeling
The mathematical model serves for solving water pressures 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 parameters.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 m 2 .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 k t z = = j 08 .; all these coefficients are dependent on z, which is axial -vertical -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 publications 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 prevalently inertial in this case.
To ensure an appropriate material for the stoppers, a special 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 engineering 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 depends on double electric layers between distinctive clay particles.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 consequent increase in the impermeability of the stopper can be attained by flocculation admixtures, which speed up the process 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 flocculation 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 average size of the grains, Fig. 4.This relationship determines the design of the construction technology of the stoppers.5 depicts the assumed time dependency on the porosity.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 flocculation, 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 significantly in the final stage, e.g., Young's modulus E dropped by 50 to 60 percent and toughness by 30 to 40 percent.

Desai's model
Desai's model is introduced because of the different volume 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, A s .Force F f is transmitted through the fluid area A f .Then the force equilibrium gives: Dividing the last relation over the nominal area, A, of the RVE, after rearrangement we obtain: This denotation from Desai's composite theory, [14], substitutes as: s s s The latter denotation is advantageous from certain points of view.Since the volume fractions sometimes change basi-  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 (linearly, 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: where h is the distance from the upper level of the stopper downwards, g is volume weight and s z,0 f is the pressure on the upper boundary of the insulator.

Some results and discussion
We first introduce the very important distribution of coefficients of filtration k (m/s) in the vertical direction depending on moisture w (%), which was derived from large-scale models.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 considered as: S w = 092 ., S a = 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 densities of the constituencies are given by: r w =10 kN/m 3 , r a = 00214 .kN/m 3 , r s =170 kN/m 3 .
The equations that have to be solved are listed below, with the coefficients mentioned above: The .
The second system of equations is: The distribution of water and air pressures, and the stresses from Lame's equations are computed by finite element method.At certain points, strength in the sense of the Mohr-Coulomb hypotheses (ideally elasto-plastic law) is attained 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 water 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.

Conclusions
Based on the results of experiments and numerical analysis, 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 situation 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.• 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 sections 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 negligible.

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 has to be changed by introducing a transformation x r = cos j, y r = sin j. (34)

It obviously holds
Inverse formulation to (21) gives:   (39) So far, an isotropic harmonic operator has been considered.For an anisotropic harmonic operator a similar approach can be applied, but natural transformation of the coefficients of anisotropic permeability In the case of axisymmetry with respect to the axis of the shaft we obtain in cylindrical coordinates in the sense of distributions:

Fig. 1 :
Fig. 1: Cross-section of the shafts and insulators -original state after abandoning the mine -geological conditions in the overburden of the mine 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 .

)
In (16) the term on the right hand side represents the inflow 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 increase in the volumetric strain of the soil skeleton & e vol by an increase in the density & r 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:

Fig. 2 :
Fig. 2: Alternative setting of more than one stopper

Fig.
Fig.5depicts the assumed time dependency on the porosity.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 flocculation, 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 of fluid and solid.

Fig. 7 :
Fig. 7: Relation coefficient of filtration k and moisture w , as the transformation itself covers the coefficient of rotational anisotropy.Substituting the above formulas to the anisotropic harmonic operator in 2D yieldsIn the isotropic case, k k k k k x y r = = = = j, depending on the position and the time.k j and k r are the permeability coefficients in the hoop and radial directions, respectively.

Table 1 :
Material properties of the spoil soil first equations read:

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