Mathematical Analysis of Transport of Contaminants Through Unsaturated Porous Media in Applied Soil Column Experiments

 

Mathematical Analysis of Transport of Contaminants Through Unsaturated Porous Media in Applied Soil Column Experiments

Rekha J1

Department of Mathematics .

Cambridge institute of Technology, Bengaluru ,India

Dr Suma S P3

Department of Mathematics .

Cambridge institute of Technology,Bangaluru India

Dr. Ramesh.T2

Department of Mathematics .

Cambridge institute of Technology, Bengaluru India

Dr.S R .Sudheendra4

Department of Mathematics .

Presidency University Bengaluru India

Abstract The main objective of the study is to provide mathematical model for better understanding of transport of pollutant through unsaturated and saturated porous media. A mathematical model is an important tool and can play a crucial role in understanding the mechanism of groundwater pollution problems. It is a simplified description of physical reality expressed in mathematical terms. Mathematical models that attempt to simulate atmospheric processes involved in groundwater pollution are based, in general, on the equation of mass conservation for individual pollutant species. Such models relate in one equation the effects of all the physical aspects and dynamic processes that influence the mass balance on groundwater which include transport, diffusion, removal of pollutants and loss or transformation through chemical reactions.

Keywords Integral Transforms Method, Mathematical Modeling Advection Dispersion Equation, Saturated porous Medium unsaturated porous media, heterogeneous flow.

I. INTRODUCTION

Solute transport in the unsaturated flow is determined by the phenomena of advection and hydrodynamic dispersion, both of which depend on the flow field. Hence, a general model of solute transport must include the equations for flow and solute displacement with terms appropriate for physical, chemical or biological interactions of the solute with the soil environment. The increasing demand for water for domestic, industrial and agricultural purposes is placing greater emphasis on the development of ground water resources. The exploitation of ground water resources at some parts of the country induces degradation of groundwater quality as well as the discharge of untreated effluents which add contaminants to the groundwater system.

Water is the most precious resources provided by the planet earth to humanity. Water forms an essential ingredient for the existence of all forms of life, and is, therefore considered to be the life blood of biosphere, and controlled growth of population and rapid industrialization have resulted in increased demand for water in human activities. Water is available on the earth either on the surface of the ground

[surface water] or below the ground surface [groundwater]. The importance and management of these two types of resources have been studied independently by many researchers over the past two decades. Now days researchers have felt the need to study both the systems together for better understanding of the interaction between surface water and groundwater in the unsaturated zone, which lies above the water table and up to the ground surface linking the two types of resources mentioned above.

For mathematical models, the earlier works made use of Green’s theorem and stochastic process for solving the partial differential equations, resulting in complicated solutions. In the present work, we have made use of the integral transforms, moving coordinates and Duhamel’s theorem to solve this equation, which simplifies the solution to a simple form, considering different parameters like adsorption, non-adsorption, chemical reaction, and radioactive decay. For the same model the results are compared with the experimental data.

Modeling studies

The mathematical models consist of integrated hydrological and solute transport models and are used to simulate contaminant movement and subsequent uptake by plants.

The hydrological model is based on Richardss equation for the transport of water through soil in response to the forces of gravity and water pressure. Richards equation relates the change of soil moisture content to the hydraulic flux divergence. It assumes that water flow can be represented by Darcys law whereby the flux is related to the product of the hydraulic head gradient and a non-linear hydraulic conductivity parameter which depends on the moisture content of the soil.

q = -K (dh/dz) q = soil water flux K = hydraulic conductivity

h = hydraulic head = elevation

The change in soil moisture content is related to the divergence of the hydraulic flux

Minus any uptake by the plants and this follows the law of conservation of matter.

d/dt = dq/dz – Uw

= soil moisture content = water uptake rate

Modeling the movement of solutes through soil is carried out through the use of the Advection-dispersion equation where the rate of change of contaminant stored within a unit volume is related to the divergence of the advective flux due to bulk water motion,(Jq) and a dispersive flux (Jd) resulting from the combination of molecular diffusion and Mechanical dispersion.

d(c + S) / dt = d(Jd+Jq) / dz – Us – T

c = activity of contaminant stored in solution S = activity of soared contaminant

Us = root uptake,F = radioactive decay

The contaminant stored in the soil is partitioned between the solute and soared phases. The sorption of contaminant onto the surface of the soil is assumed to be represented by

a linear partition coefficient (Kd) according to the following expression. S = Kdc

The uptake of solute by the plant root system is represented using a linear form of Michaelis-Menton rate kinetics as described by Nye and Tinker in 1977.

Us = 2praa Us = uptake by the plant

Fate & Effects of Pollutants:

One major influence is on the abundance of films, and is strictly geometric. From dimensional considerations of geometrically similar porous media, both the specific surface (surface area per unit volume), and specific perimeter (Ps, grain perimeter length per unit crosssectional area) scale as 1. Therefore, at low water saturations, thin film flow is localized on grain perimeters, which per unit bulk area increases as grain size decreases. A functional relation between Ps and can be obtained for mono disperse packs of spherical grains by first determining the number of grains N0 intercepting a plane of unit area A0 in a porous medium. Such a plane will intercept grains with centroids located within a distance of /2 on either side of it. Therefore, grains centered within a reference bulk volume of unit crosssectional area A0 and thickness will have perimeters that project onto the unit crosssectional surface. With a porosity n, the solid phase volumetric fraction of 1 n is equal to the ratio of the solid volume per unit bulk volume. For the reference volume of A0,

Modelling of the filtration

In the modelling of the underground transport, first question is how to handle the fluid flow. In all models reviewed in

[23] and [15], and further studied in [14], we remark that the concentration of the transported contaminants is such that the flow is not affected. This means that we should study the

effective flow before looking at the mass balance equations for the transported species. All flows could be considered as incompressible.

For the multiphase flows, we can’t say much about the derivation of the microscopic models from the first principles. We refer to the review paper of Bourgeat [12] and to the calculations for the flow through a singe pore in [19]. To the contrary, homogenization of the microscopic models was studied with success in many papers. We refer to [9] and

[18] for the 2-phase flows in periodic geometries and with the permeability of same order, to the article [10] for the random case and to [11] for the 2-phase dual-porosity model.

Very frequently, solutes are transported by an unsaturated flow. Single phase unsaturated flows are described by the Richards’ equation. Even if it seems that the equations describing the 2-phase flow of water and air are more realistic, the Richards’s equation is widely used and in the hydro chemical computer codes one should solve it in a heterogeneous medium.

Hence in order to have an effective equation for the velocity and the pressure, we study the homogenization of the Richards equation. The Richards equation in the dual porosity case reads as follows (see e.g. [7])

Where d is a small positive number, Pe, pe are the unknown (scaled) pressures, J is the water content and k is the hydraulic conductivity. Let us note that we have two functions J and k. Each J is a C1 strictly monotone increasing positive function on , equal to ds for s > 0. Each k is a C1 increasing function defined on [0,1]. They are positive and strictly monotone increasing for J > q0. They could be zero on [0, q0 ). Let V = { v ÃŽ H1 (W) | v = 0 on G1 } .

Proposition 1. Let PD Î H1 (0,T; L¥ (W)) Ç L2 (0,T; H1 (W)),

PD > J-1 (q0 ) +d, d > 0. Let g Î L2 (G2 ×(0,T)). Then there

exists at least one solution e = Pe + Pe ÃŽ (PD + L2

(0,T; V)) Ç L¥ (W×(0,T)) , e > J-1 (q0 ) + d and

+ ÃŽ L2 (0,T; V’) satisfying

and satisfying the initial condition (10). Furthermore, it satisfies the a priori estimates

where constants are independent of e.

Proof. Proof of the proposition follows the lines of the article [5]. It should be noted that the conditions on the data guarantee the non-degeneracy of the coefficients.

Next we use the results from [1] presented in the introduction and extend Pe to W. The extension PePe satisfies

Lemma 2. We have

Using the a priori estimates derived in proposition 1 and lemma 2 and the concept of two scale convergence as in [3], we get the following compactness results

Proposition 3. There exists a subsequence such that

Proof. The proof is a direct consequence of the a priori estimates and of theorem 1.2 and proposition 1.4 in [3].

Again as in [3] we plug into (11) a test function of the form j(x,t) + ej1 (x, x/e, t) + Y(x, x/e, t), where j ÃŽ (QT), j1 ÃŽ (QT ; (Y)) and Y ÃŽ (QT ; (Y)) with Y = 0

for y ÃŽ Wf. Passing to the two-scale limit yields:

where Qm is the two-scale limit of the (sub)sequence of Jm (pe) . Let xj, òYf xj = 0, be the 1-periodic solution to

and X the tensor whose (i, j) component is the average of ¶xj

/¶yi . Then P1 = åj xj (y) ( ¶ P / ¶ xj-dj3 ) and the choice Y = 0 gives

We note that by Lemma 4.2 from [11], ¶t òYmQm dy Î L2

(W×(0,T)).

The real difficulty is in identifying Qm. Let us note that pe doesn’t converge strongly and oscillations are likely to persist. We proceed by using the the method of periodic modulation as in [11]. It was introduced in [6] and we repeat the definition:

Definition 4. For a given e > 0, we define a dilation operator De mapping measurable functions on × (0,T) to measurable functions on W × Ym × (0,T) by

where ce (x) denote the lattice translation point of the e – cell domain containing x. We extend De u from Ym to Èk (Ym + k) periodically.

It is easy to see that

Lemma 5 For “u ÃŽ L2 (0,T ; W1,q ( )) we have

Furthermore, let j, y ÃŽ L2 (0,T ; ). Then

Proposition 6 ([11]). Let {ue} be a uniformly bounded sequence in × (0,T)) which satisfy the conditions

Deue u0 weakly in L2(QT (Ym)) and ue ® u* Î L2(QT ; (Y )) in the two-scale sense.

Then we have

Proposition 7. With the above notations we have

Let e > 0 and let De pe be defined by (25). Then De pe = e

satisfies in L2 (0,T;H-1 (Ym)) the equation

We consider the problem:

For fixed P Î L2 (0,T ; H1 (W)) Ç L¥ (QT) ,J-1 (q0 ) +d < P (a.e)

on QT and P0 Î L¥ (W) ,J-1 (q0 ) + d < PD (·, 0) (a.e) on W, the theory from [5] and [21] gives existence of a unique solution p* for (34)-(35).

Our final step is to compare the problems (34)-(35) and (32)- (33). We have the following result:

Proposition 8. Let De se be defined by (25) and let p* be a solution for (34)-(35). Then we have

where p is defined by (19).

Proof. We choose the test function by a change of pivot. Therefore, we introduce for (a.e.) (x, t) ÃŽ QT a function we ÃŽ (Ym) by

Obviously, we Î H1 (QT×Ym), ||Ñywe (x,y,0)

|| < Ce and

Now we introduce the function ze by ze = Jm ( De pe ) -Jm ( p* ) and Km (s) by Km (s) = km(Jm (h)) dh. Obviously, ze satisfies the vibrational equation

and the boundary condition

ze = Jm (De Pe) – Jm(P) in L2(0, T; H1/2(¶Ym)) for (a.e.) x ÃŽ W.

In order to estimate the L2 – norm of ze we choose Y = we as the test function. We have

Using the uniform bound for Ñywe in L2(QT × Ym) we find out that the right hand side is bounded by C e.

Furthermore

Therefore, we write (39) in the form

It remains to find an estimate from below for the diffusion term. We have

Therefore

and consequently,

{Km(De pe) – Km (p*)}{Jm(De pe) – Jm(p*)} ® 0 (a.e.) in QT

×Ym.

Hence, after passing to a subsequence, Depe – p* ® 0 (a.e.) in

QT × Wm and Depe – p* ® 0 (a.e.) in Lq ( QT × Ym) for “q ÃŽ [ 1,

+Â¥[ by Lebesgue’s dominated convergence theorem.

Therefore the sequence { Pe, pe } converge to the unique1 solution { P, p } of the problem

3 MODELLING TRANSPORT, DIFFUSION AND ADSORPTION OF CHEMICAL SPECIES

We are going to assume that certain chemical species are dissolved in the fluid. Their concentration is sufficiently small and the fluid flow is not affected. For simplicity of notation we write the equations for just one specie; the generalization to m substances is straightforward. The conservation of mass law is

where ve = kf(Jf (Pe)) (ÑPe – e3 ) is the filtration velocity, ve ue is the convection term, -DÑue is the dispersion and diffusion term and f(ue) represents the chemical reactions.

The simplest way to include adsorption effects is to introduce the concentration Ue of the substance absorbed on the pellet surface. It is reasonable to suppose that the actual adsorption rate be depends on ue and Ue. Thus we have

where b is a Holder continuous function of u and U. In order to complete the model we add an equation describing the dynamics of the adsorbed concentration Ue which reads

Using the results from [17] we find out that { ue, Ue } converges in 2-scales to the solution { u, U } of the following problem

We suppose

g1 Î L2(G2 × (0, T)) and Ug, ug Î H1 (0, T; L¥(W)) Ç L2(0, T; H1 (W)).

Well-known examples of nonlinear adsorption rates are the Freundlich isotherm b(u,U) = g( up -U), g > 0, 0 < p < 1 and the Langmuir isotherm b(u, U) = g(au / (1+bu) – U), a, b, g >

0. The next step would be to add the surface diffusion as in

[16] or to consider the nonlinear transmission to the pellet, as in [17]. Getting more complicated models from [23] would be subject of our subsequent publications.

4 NUMERICAL SIMULATIONS

We would like to point out the effects of chemical reactions. Consequently, without losing generality, we consider only a saturated flow in W = (0, 2p)2 and the periodic boundary conditions for the concentrations. To obtain a precise approximation to the saturations, we approximate the spatial derivatives using the pseudo-spectral method. More precisely, we write u and U in terms of the corresponding spectral Fourier basis. The unknown functions in the nonlinear terms are evaluated in the collocation points, coming from Fourier’s grid of the size 256 ×256. Then using the inverse FFT, the nonlinear terms are written in the spectral basis.

In the time discretization, we integrate exactly the linear terms: (a) the diffusion-convection term in (50) and (b) the linear term containing U in the adsorption rate b. For the nonlinear terms we use a third order explicit Adams- Bashforth scheme for f(u) in (50) and an implicit third order Adams-Bash forth scheme for b(·, u, ·) in (54). After scaling the equations with the characteristic time Tc = 2×104, we take the time step 10-3.

This numerical approach was developed in 2D turbulence simulation by Robert and Rosier (see [22]) and the method is very robust for the non-linear coupled transport and diffusion. We present two numerical experiments. They are with same data, but with different isotherms. The data are given by the following table:

f(u) is taken to be linear. The results are at the figures which follow. They show important effects of the transport on the evolution of the concentration profiles. Also diffusive effects depend very much on the choice of the isotherms, as visible on the figures.

After these calculations, which correspond to the simplest homogenized model, we are planning to develop a numerical method for the unsaturated flow, to add the surface diffusion in the modelling of the adsorption and to include the effects of high Péclet number. The basic challenge is to develop fast solvers for dual-porosity models and for more details we refer to the recent article [2].

Waters in Mono Basin, California (Neumann and Dreiss, 1995). A follow-up study was presented to remeasure and reevaluate the 36Clcontent at several previously studied sites in the Great Basin (Phillips et al., 1995). Soil-water chemistry

and isotopic data were presented from three deep vadose zone boreholes at the Nevada Test Site to quantify soil-water flux and its relationship to climate (Tyler et at., 1996). Tracing techniques and core data were used to study the temporal and partial variations of the soil-water flux and recharge during the last 120000 years. Natural variations of stable isotopes were used to study infiltration in clay soils and groundwater recharge in the small drainage basin of Barogo in West Africa (Mathieu and Bariac, 1996). The isotope analyses led to the characterization of two separate infiltration processes as well as two recharge processes.

Stochastic Modeling and Analysis.

The spatial-temporal averaging method was used for modeling transport in porous media with a nonhomogeneous distribution of elementary domains in the spatial-temporal space (He and Sykes, 1996). Several averaging theorems and corollaries about the averages of spatial and temporal derivatives were presented and rigorously proved .The real- space expressions for second-order mean velocity and velocity covariance in 3-D, statistically anisotropic media were derived, and the effect of higher-order flow correction terms on First- and second-order reliability methods were applied as alternatives to Monte Carlo simulation in the probabilistic analysis of groundwater contaminant transport and remediation (Hamed et al., 1996). A 2-D FE model was interfaced with a reliability analysis program that provided the probability that a contaminant exceeds a target level at a well. Covariances of velocity, head, and log (hydraulic conductivity) under quasi steady flow were developed, and their effect on adventive transport was examined (Zhang and Neuman, 1996a). Results showed that temporal fluctuations in the direction of the mean velocity cause longitudinal dispersion to decrease and transverse dispersion to increase. Probability distributions of travel times were developed for 1- D (vertical) convective solute transport in field-scale soils under unsteady and non-uniform flows .Approximate ensemble probability distribution functions of conservative solute travel time were derived directly from the convective transport stochastic partial differential equation. A higher- order theory was presented for steady-state mean uniform saturated flow and nonreactive solute transport in randomly heterogeneous porous media (Hsu et al., 1996). Results show significant improvements compared with lower-order theories. Dagan’s solution for transport in heterogeneous porous formations was extended to the case of finite Pe (Fiori, 1996).

It was suggested that pore-scale dispersion matters only with reference to transverse spreading, and that Dagan’s solution, valid for Pe = 00, was an adequate approximation in a wide range of finite Pe. General equations for the first two moments of the particle trajectory were derived for the problem of transport of inert solutes in random non stationary velocity fields . It was shown that the equations can be solved analytically for the case of quasi-unidirectional mean flows, and the theory was applied to the case of transport in media displaying a linear trend in the mean log (conductivity). The general expressions for the time-dependent ensemble averages of the second spatial moments (A) and the effective dispersivities ‘Y were evaluated to study the effect of initial plume size on (A) and ‘Y in 3-D heterogeneous isotropic

aquifers (Zhang, Y.K., et at., 1996). Results confirm that (A) and ‘Y approach their respective erotic limits as the size of the source increases and that transverse lengths of a source were more important than longitudinal length for the erotic condition to be met. The exact and approximate deterministic partial differential equations of the time-space evolution of mean, two-point moment, and n-point moment solute concentration for conservative solute transport by steady and unsteady groundwater flow in a heterogeneous aquifer were developed in the real space time domain. The derivations were performed by means of the cumulant expansion method, combined with the calculus for the time-ordered exponential and with the calculus for the Lie operator. A fully nonlinear analysis of the integral-differential equation for the spatial second moments X of the ensemble mean concentration in a heterogeneous aquifer was carried out by solving the nonlinear integral-differential equation for X numerically and iteratively (Zhang and Chi, 1995). The effects of log (hydraulic conductivity) variance, Pe, and anisotropy then were investigated.

.

Hydraulic properties of adsorbed water films in unsaturated porous media

Hydraulic properties of adsorbed water films in unsaturated porous media, Volume: 45, Issue: 6, First published: 12 June 2009, DOI: (10.1029/2009WR007734)

Hydraulic properties of adsorbed water films in unsaturated porous media

Hydraulic properties of adsorbed water films in unsaturated porous media, Volume: 45, Issue: 6, First published: 12 June 2009, DOI: (10.1029/2009WR007734)

Hydraulic properties of adsorbed water films in unsaturated porous media

Hydraulic properties of adsorbed water films in unsaturated porous media, Volume: 45, Issue: 6, First published: 12 June 2009, DOI: (10.1029/2009WR007734)

CONCLUSION:

A discussion was presented regarding a paper on inappropriate use of analytical fate and transport models to estimate the A stochastic modeling approach and Monte Carlo simulation were used to characterize the uncertainty in migration pathways in geological media surrounding a hypothetical buried water repository in Bedfordshire, U.K. (Mackay, R., et at., 1996). The coupled effects of porous media heterogeneity and flow unsteadiness on transport were investigated using a stochastic- Lagrangian approach, where log(conductivity) was modeled as a space random function and flow unsteadiness was considered deterministic, affecting the magnitude and direction of mean head gradient (Dagan et at., 1996a). By applying the model to the case of a periodic variation in the flow direction, it was shown that unsteadiness has an insignificant effect on the longitudinal spread, whereas it strongly affects the transverse spread. An experimental velocimetry technique based on particle image analysis was used to perform a Lagrangian description of nonreactive pollutant particle motion in a 3-D saturated porous medium (Cenedese and Viotti, 1996). Statistical analysis of the experimental data allowed for estimation of velocity and displacement probability density functions, velocity component correlation functions, Langrangian integral scales, and mechanical dispersion coefficient tensor components. A recently developed exact Eulerian-Lagrangian theory of advective transport in space-time random velocity fields was extended to account for anisotropic local dispersion (Zhang and Neuman, 1996b).

REFERENCES:

1. Atul Kumar Verma, Murthy Bhallamudi, S. and V. Eswaran (2000): Overlapping control volume method for solute transport, J. of Hydrologic Engg., : 308-316.
2. Ahuja, L.R., De Cousey, D.G., Barnes, B.B, and K.W. Rajas (1993): Characteristics of macropore transport studied with ARS root zone water quality model. ASAE. 36-2, March-April, pp. 369-380.
3. Al Niami, A.N.S., and K.R. Rushton (1979): Dispersion in stratified porous media: Analytical solution. Water Resour. Res., 15 (5): 1044- 1048.
4. Barry, D.A., Parlange, J.Y., Sander, G.C., and M. Sivaplan (1993): A class of exact solutions for Richards equation. J. of Hydrology. 142: 29-46.
5. Basak, P. and V.V.N. Murthy (1978a): Pollution of groundwater through nonlinear diffusion. J. of Hydrology. 38: 243-247.
6. Basak, P., and V.V.N. Murthy (1981): Groundwater quality improvement through nonlinear diffusion. J. of Hydrology. 53: 151- 159.
7. Bear, J (1972): Dynamics of Fluids in Porous Media. American Elsevier, New York.
8. Bear, J., and A. Verruijt (1990): Modeling Groundwater Flow and Pollution. D Reidel Publishing Co., Tokyo.
9. Bitton, G., Davidson, J.M., and S.R. Farrah (1979): On the value of soil columns for assessing the transport pattern of viruses through soils A critical outlook. Water, Air and Soil Pollution. No. 12, pp. 449-457.
10. Brandt, A., Bresler, E., Diner, N., Ben-Asher, J., Heller, J., and D. Goldberg (1971): Infiltration from a trickle source: (i) Mathematical Models. Soil Sci. Soc. Am. Proc. 35: 675-682.
11. Carslaw H.S. & Jeager J.C. (1947), Conduction of Heat in solids, Oxford University Press, P 386.
12. Fletcher (1991), Computational techniques for fluid dynamics, 2nd edition, Springer Verlag, Berlin.
13. Shamir, U.Y., and D.R.V. Harleman (1967) : Numerical solution for dispersion in porous medium, Water Resour. Res. 3 (2): 557-581.
14. Thomas, G.W., and A.R. Swoboda (1970): Anion exclusion on chloride movement in soil. Soil Sci., Vol. 110, pp. 163-166.
15. Tisdale, S.L., Nelson, W.L., Bealon, J.D., and J.L. Havlin (1997): Soil Fertility and Fertilizers. Prantice Hall of India Pvt Ltd., New Delhi.
16. van Dam, J.C., & R.A. Feddes (2000) : Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation. J. Hydrology. 233:72-85.
17. Zheng, C. and Bennett. G.D. 1995: Applied Contaminant Transport modeling, Theory and Practice. Van Nostrans-Reinhold. NY. U.S.A.