A molecular approach to the foundations of the theory of solute transport in porous media

Journal o f Hydrology, 50 (1981) 355--370

355

Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands [2] A MOLECULAR APPROACH TO THE FOUNDATIONS OF THE TIIEORY OF SOLUTE TRANSPORT IN POROUS MEDIA I. C o n s e r v a t i v e S o l u t e s i n H o m o g e n e o u s I s o t r o p i c S a t u r a t e d M e d i a

VIJAY K. GUPTA 1 , R.N. B H A T T A C H A R Y A z and GARRISON SPOSITO 3'*

1 Department o f Civil Engineering, University o f Mississippi, University, MS 38677 (U.S.A.) 2 Department o f Mathematics, University o f Arizona, Tucson, A S 85721 3 Department o f Soil and Environmental Sciences, University o f California, Riverside, CA 92521 (U.S.A.) (Received February 15, 1979; revised and accepted May 22, 1980)

ABSTRACT Gupta, V.K., Bhattacharya, R.N. and Sposito, G., 1981. A molecular approach to the foundation of the theory of solute transport in porous media, I. Conservative studies in homogeneous isotropic saturated media. J. Hydrol., 50: 355--370. A molecular model is developed for the transport of a conservative solute at low concentrations in a homogeneous isotropic water-saturated porous medium. In this model, a solute molecule undergoes contact collisions with the solid grains in the medium at successive random times; in between these collisions, the velocity of the molecule is governed by the Langevin equation; the effect of the collisions with the solid grains is to scatter a molecule in a random direction. This molecular model is employed to derive rigorously a parabolic differential equation for the solute concentration at the macroscopic level. In the absence of solute convection, the coefficient of molecular diffusion in a porous medium is proved to be less than the coefficient of molecular diffusion in bulk solution, a finding which is in agreement with experimental observations. For non-zero convection, the assumption of isotropicity of the medium is employed to prove that the solute dispersion tensor is diagonal. For small magnitudes of the liquid velocity, the isotropicity assumption also implies that the coefficients of longitudinal and transverse dispersion are approximately parabolic functions of the liquid velocity. The expressions derived for the dispersion coefficients, when compared with their experimentally observed values, suggest that their dependence on the liquid velocity comes primarily through solute--liquid molecular collisions instead of through collisions of the solute molecules with the grains of the solid phase. The solute convective velocity is shown to be less than or equal to the liquid velocity. Precisely when the difference between the two velocities will be significant remains to be established.

* The ordering of authors' names in this cooperative investigation is not significant.

0022-1694/81/0000--0000/$02.50

© 1981 Elsevier Scientific Publishing Company

356 INTRODUCTION

The isothermal transport of a conservative solute through a water-saturated homogeneous porous medium, at the macroscopic level, is commonly taken to be governed by the partial differential equation: 3

~2 C

~C _ ~_, Dij ~)t i,j=l ~)xiSxj

3

O(UiC)

,V i=~

3xi

t :~ ~

where C ( x , t ) is the solute concentration (M L-j ) expressed as mass per unit volume of medium; D~j ( L 2 T -1 ) is an element of the solute dispersion tensor,/9; and U is the mean convective velocity (L T -1 ) of the fluid phase in the porous medium. Eq. 1 is obtained after combining the equation of continuity: OC

~-~-+ ~v'(Jdis

+

_Jco.) = 0

(zi

with assumed empirical expressions for the solute mass flux density vectors given by: Jd~ = - - D ' V C

and

~.]con --- UC

(3a, bl

At present, the principal justification for the use of eq. 1 in the study of dispersion p h e n o m e n a comes directly from experimentation (see, e.g., Bear, 1972). Although at times discrepancies between the solutions of eq. 1 and experimentally observed concentrations have been reported in the literature and tentative modifications of eq. 1 to account for these discrepmacies have been proposed (see, e.g., Fried and Combarnous, 1971), eq. 1 remains the, most widely used equation in macroscopic solute transport studies. During recent years, the derivation of eq. 1 from fundamental physical considerations has become an important research area (Beran, 1968; Fried and Combarnous, 1971; Bear, 1972; Scheidegger, 1974). The expectation motivating these fundamental investigations, as well as the present one, is that a mathematically rigorous, physically meaningful derivation of eq. 1, or a closely related differential equation, will lead to a deeper understanding of the limits of applicability, in the field, of the currently accepted macroscopic solute transport theory. In this paper, a derivation of an equation very similar to eq. i is presented that is based on a dynamical model of solute transport at the molecular level. The molecular approach to the foundations of macroscopic solute transp o r t equations has been discussed and elaborated on by Beran (1968), Scheidegger (1974) and Todorovic (1975, and the references cited therein). The theories and model calculations of these investigators, as well as those of other workers, have been reviewed critically very recently by Sposito et al. (1979). In this review, it was shown that a principal difficulty with the foundational theories that have been developed previously at the molecular level is their c o m m o n lack of a dynamical basis. These theories often state,

357 correctly, the basic need to describe the stochastic velocity process of a solute molecule moving through a porous medium and to compute, ultimately, the asymptotic expression for the probability density function (p.d.f.) of the corresponding position process of the molecule. But all have stopped short of establishing a realistic dynamical model of solute movement that could lead in a rigorous fashion to equations for the velocity and position processes and then on to a derivation of eq. 1 or a related equation. The theory developed in the present paper commences with dynamicallybased expressions for the change in m o m e n t u m of a conservative solute as caused by its interactions with water molecules and by its collisions with the grains of the solid phase in a water-saturated porous medium. The effective mass of the solute molecule is assumed to be large in comparison with the mass of the water molecules. The solute molecule may be electrically neutral or charged, but regardless of this, it is assumed to be only in weak interaction with other solute molecules and consequently, the transport of immiscible substances or of miscible substances at high concentrations cannot be described by the present model. On the other hand, the model should apply to the solute transport phenomena encountered in soils and aquifers whenever the concentration of solute is low enough to make the solute dispersion tensor independent of the solute concentration, as it is in eq. 1. In the following section, the specification of the equations governing the change in mom e n t u m of a solute molecule will be described in physical terms. The solution of these equations completely determines the stochastic velocity process {V(t); t > 0} of a solute molecule. The principal outcome of this construction is the fact that, in an appropriate asymptotic sense, the position process {X(t); t > 0}, as a time integral of the constructed velocity process, is asymptotically Markovian (a Brownian motion) with a mean drift vector, ~(U), and the dispersion matrix, D(U). Because this result is asymptotic in nature it is valid at the time and space intervals over which a solute molecule undergoes a very large number of collisions with both the liquid molecules and the grains of the solid phase, i.e. at the macroscopic space and time scales. The Markovian nature of the position process, along with the assumption of weak interaction among solute molecules, then may be used to show rigorously that the macroscopic solute concentration, C(x,t), is a solution of a parabolic partial differential equation. This equation differs from eq. 1 in t h a t the mean convective velocity of the solute, ~(U), which appears in it is n o t necessarily equal to U, the mean fluid velocity. In a separate section, a detailed examination of the expression for the mean solute velocity, ~(U), shows that for some porous media, (e.g., unconsolidated media, ~(U) m a y be approximated well by U. For this case, the derived parabolic equation governing the solute concentration reduces to eq. 1. Although ]~(U) I ~ t U ] if U ¢ 0, it is not clear at present how significant this difference between the two mean velocities will be for any given porous medium.

358 A precise understanding of this important aspect will require careful experimentation with various porous media as well as some further theoretical analyses. When U = 0, the coefficient of dispersion matrix D0(U) reduces to the coefficient of molecular diffusion matrix Do. The derived expression for the coefficient of molecular diffusion in porous media is show~ to be less than the coefficient of molecular diffusion in bulk solutions. This f!nding is in agreement with experimental observations. I t is also shown that D(U) is diagonal. For small values of U, the elements of D(U) can be approximated by quadratic functions of U. Furthermore, the derived expressions for the diagonal elements, i.e. the coefficients of longitudinal and transverse dispersion, when compared with their experimentally observed values suggest t h a t their dependence on the fluid convective velocity U comes primarily through collisions of the solute molecules with the water molecules instead of their collisions with the solid phase. The paper is concluded with some observations on the scope of further research on the foundations of solute transport in porous media.

SPECIFICATION OF THE DYNAMICS GOVERNING THE VELOCITY OF A SOLUTE MOLECULE Consider a macroscopically homogeneous isotropic water-saturated porous medium. By homogeneity it is meant specifically that the dynamics of transport are invariant under all translations of the coordinate axes fixed in the solid phase. In order to understand the term " i s o t r o p y " , imagine two rectangular Cartesian coordinate systems, X,Y,Z and X',Y',Z', having the same origin. Then isotropy means that the expression for the general dynamical equations (including their functional forms and the parameters involved therein) with respect to the two coordinate systems is identical (excluding terms related to the external force of gravity). Let the initial time of observation be t -- 0, and let C o (x) be the solute concentration in the porous medium at this initial time. It is assumed that the mean velocity of the liquid, U is steady and uniform; t h a t the solute is conservative (i.e. nonadsorbing and nonprecipitating) and is in low concentration; that the transport is under isothermal conditions, and that, for all times t ~ 0, the total solute mass in the porous medium remains fixed. A typical experiment depicting this situation may be conceived as follows. Consider a vertical column filled with sand grains of nearly uniform diameters, whose arrangement satisfies the conditions of macroscopic homogeneity and isotropicity. A steady, uniform flow of solution through the column of sand can be produced by installing a tank above the column and maintaining a constant solution depth in the tank. As a conservative solute injected at some point in the medium in low concentration, one can consider the radioisotope 131I in the form of Nal (see, e.g., Todorovic, 1970). The specifications of the solute dynamics at the molecular level now will

359 be considered in detail. The assumption that the solute is in low concentration implies that the solute molecules are in weak interaction among themselves. Consequently, the change in m o m e n t u m of a solute molecule can be derived solely from its interactions with water molecules and with the solid phase. The statistical implication of this fact, in turn, is that the velocities of solute molecules may be treated as being stochastically independent of one another. In the subsequent development, it is assumed that the molecular interactions between the solute molecules and the grains of the solid phase are essentially zero-ranged and henceforth, for this case, the terms "molecular interactions" and "molecular collisions" will be used interchangeably. The assumption of zero-range forces is employed c o m m o n l y in the kinetic theory of gases (Present, 1958). We begin our discussion with the Langevin equation, which at an "appropriate time scale" (a concept to be discussed later) describes the rate of change of m o m e n t u m of a solute molecule moving in the liquid phase (Chandrasekhar, 1943; Rice and Gray, 1965). The Langevin equation provides the stochastic dynamical basis, at the molecular level, of diffusion in liquids, i.e. of Fick's first law. In the present context, although the solid phase introduces a significant change in the dynamics, the Langevin equation still plays an important role in the construction of a stochastic dynamical model of solute movement in a porous medium. For this reason, we shall first discuss, only heuristically, what the Langevin equation represents and how it comes about. Consider a typical solute molecule moving in the liquid phase. The force causing this m o t i o n results from incessant collisions (i.e. short-ranged interactions) between the solute molecule and the liquid molecules. Typically, the intervals between'collisions for very large molecules, e.g., colloids, are of the order of magnitude of 10 -2~ s and for typical hydrated ions are of the order of magnitude of 10 -14 s. On this time scale, the forces experienced by the solute molecule exhibit strong statistical dependence, because of the many-body nature of the interactions of the solute molecule with the water molecules. However, as the time evolves, the dependences between forces which are applied sufficiently apart in time (in relation to the collision times) become weak (see, e.g., Rice and Gray, 1965). Now, if one considers the net change in the m o m e n t u m of a solute molecule, n o t as produced by individual collisions but instead as produced by the sum of a large number of collisions with the water molecules, one may expect the forces causing the net change in m o m e n t u m to be essentially independent in time. As a result of this condition, one can also expect the velocity process of the solute molecule to be governed by a Markov process. Since the net m o m e n t u m change is a sum of the individual changes in m o m e n t u m , one can expect the net mom e n t u m change to be Gaussian; this result essentially follows from the central limit theorem (see, e.g., Bhattacharya, 1978). Indeed, according to the Langevin equation (applicable over a time scale which is large in comparison with the collision times), the velocity is governed by a Gaussian, Markovian

360 stochastic process. A detailed discussion of the mathematical structure of the Langevin equation and the time scales over which it describesphysically the solute velocity will be given later (see also Rice and Gray, 1965, ch. 4). Now, for a typical solute molecule moving through a water-saturated porous medium, it is reasonable to assume that the solute molecule collides much more frequently with the water molecules than with the solid phase~ This assumption will enable one to determine separately the change in mom e n t u m of a solute molecule caused by these two types of collisions and then to combine them together to determine the structure of the velocity process. Specifically, let 0 ~ To <: T1 • • • <: ~, •. • be a sequence of random times at which a solute molecule undergoes collisions with the solid phase,, the latter now being regarded as an immobile assembly of infinitely hard grains. Three postulates will be made. The first is: (A1) The time intervals, (Ti+I -- Ti), i ~ O, form a sequence of independent and identically distributed (i.i.d.) random variables, with mean p(U). This assumption is reasonable in view of the assumed homogeneity of the medium and the invariance of the dynamics with respect to time displacements (e.g., the convective fluid velocity is steady and the gravitational acceleration, g~ is a constant). (A2) Within each interval (Ti,Ti÷,), i ~ 0, the solute molecule suffers a very large number of collisions with water molecules. Its change in momentum, mAV(t), is governed by a version of the Langevin equation, given as: mAV(t) = -- m~(U) [V(t) - U] At + m g A t + mo(U)Zk/~(t)

t4~

In eq. 4, m is the effective mass of a solute molecule, ~(U) is the parameter which scales the average interaction per unit mass between the solute molecule and the liquid molecules, g is the gravitational acceleration, o(U) is a 3 x 3 non-singular symmetric matrix, and {Bft); t > 0} is a three-dimensional standard Wiener process [Brownian motion process; see, e.g., Bhattacharya et al. (1976)]. It was mentioned earlier t h a t eq. 4 has been employed in physics to provide a stochastic dynamical basis for the phenomenon of molecular diffusion in liquids (see, e.g., Uhlenbeck and Ornstein, 1930; Chandrasekhar, 1943; Nelson, 1967). The Langevin equation has been derived by Mazur and Oppenheim (1970), starting f r o m Hamilton's equations of motion and the Gibbs distribution on the phase space o f a liquid, under the assumption t h a t certain molecular correlations ate short-lived in time and that the effective mass of the solute molecule is large compared to the mass of a liquid molecule. In the c o n t e x t of solute transport in a porous medium, this mass criterion will be m e t readily for large molecules, e.g., organic pesticides. Since the effective mass of a hydrated, monatomic ion includes the mass of the water molecules in the hydration shell around it (Chu and Sposito, 1980), as a first approximation, it may be reasonable to assume that the velocity of such a solute molecule is also governed by a version of the Langevin equation.

361 In eq. 4, the first two terms on the right-hand side on the solute molecule and the last term gives the precisely, mAV(t), the net change in m o m e n t u m in At), is approximately Gaussian (conditionally given mean: E [ m A V ( t ) [ V(t) = v]

-~ [-- mt~(.U)(v -- U) + rag] A t

describe the mean force fluctuating force. More the time interval (t, t + V(s) for 0 ~ s ~ t) with (4a)

and the covariance matrix: E [ m A V i ( t ) m A V j ( t ) [ V(t) = v] ~" [a* ( U ) ' a ( U ) ] i j A t ; V(t) = [ V l ( t ) , V : ( t ) , V 3 ( t ) ] ,

i,j = 1 , 2 , 3

(4b)

where a*(U) is the transpose of the matrix a(U), and E [ - ] denotes the mathematical expection. The parameters fi, a in eq. 4, as well as the parameter p (recall A1) may depend only on macroscopic properties of the liquid and of the medium, e.g., porosity, viscosity, temperature, U (the macroscopic velocity). At this stage we are not postulating that the parameters necessarily depend on U, but merely leaving this possibility open. In the sequel we shall, however, demonstrate that this dependence exists, it governs the magnitude of the mean convective velocity of solute and lies at the root of what is commonly referred to as dispersion in a porous medium (see the section on the structure of the transport coefficients). In the ensuing text, the possible dependence of fi and /a on U usually will not be displayed for notational convenience. The time scale over which the Langevin equation applies may be understood as follows. Let ~ be some time interval large enough (in comparison with the solute--liquid, molecular collision times) such that the stochastic forces become essentially uncorrelated over this interval, but such that ~ is smaller than fi-i. Then eq. 4 is applicable for all times t ~ ~?. A detailed discussion of this important aspect may be found in Rice and Gray (1965, ch. 4). In order to complete the specification of the dynamics, it is necessary to describe the stochastic dynamics governing the collisions of a solute molecule with the solid phase. This leads to the third postulate: (A3) The velocity, V(Ti+), of a solute molecule immediately after a collision with a solid grain at some time, Ti, is given by: V(Ti +) = p I V(Ti) lei

(5a)

where I" I denotes the magnitude of a vector, p is a constant such that 0 ~ p 1, and: V(Ti + )

= lim V ( T i + s ) s40

(5b)

Here {ei } is a sequence of independent random vectors, each having the uniform distribution on the unit sphere { Iv I = 1}.

362 In eq. 5, the vector I V(vi) l ei gives a random direction to the velocity ol ' the solute molecule w i t h o u t changing the magnitude. Stated another way, the velocity vector of the solute molecule is pointed in a random direction after the collision with the solid phase. This type of collision is known as diffusive, as opposed to specular (see, e.g., Present, 1958, p. 56). Furthermore, the solute molecule is not likely to be trapped in the solid phase, since the solute has been assumed to be conservative. The parameter p in eq. 5 accounts for a loss in the magnitude of the velocity in case a fraction of the kinetic energy of the molecule is dissipated as a result of a non-elastic eolli sion. Eq. 5 follows from solving the equations of m o m e n t u m and energy conservation governing the collision of a moving particle with a large fixed particle. The isotropy of the medium forms the physical basis for the a s s u m p tion that e~,i > 1, has a uniform distribution on the surface of the unit sphere. Note that:

E[ei]

= O,

i ~1

~6!

Eqs. 4--6 can be employed together to determine the stochastic velocity process, {V(t); t > 0}, and the position process, an integral of the velocity process. In the following section we show how this construction leads to the derivation of a macroscopic, parabolic differential equation governing solute transport, including expressions for the transport coefficients.

DERIVATION OF THE SOLUTE TRANSPORT EQUATION In this section the structure of the position process {X(t); t > 0} will be investigated and an appropriate macroscopic differential equation governing the solute concentration in a homogeneous porous medium will be derived. The position process is defined by: l

X(t) = x0 + j V(s)ds 0

In eq. 7, x0 denotes the non-random initial position of the solute molecule. The entire development to be given in this section is based on the limil theorem proved by Bhattacharya and Gupta (1979) on the basis or" the postulates (A1)--(A3). Recall t h a t / l = E [ Ti+1 -- ri ] is the mean collision time of a solute molecule with the solid phase. Write:

W = U+g/fl,

f(z) = 1 - - 1

i:(1--e-~S)Fl(ds) 0

363

~=

[

1 ;(l_e_~U~)Fl(ds) I--~o

1 W = f(3p)W

(8)

where F 1 (s) is the c o m m o n probability distribution function of the random variables (ri+ 1 -- %)/P, i > 0. Let W " W denote the matrix whose (i,j) elem e n t is WiWs, and let o ' ( U ) denote the transpose of the matrix o(U). Also write" 1 7

1

g(z) = 1 - - - - J { 2 ( 1 - - e -zs) Z

h(z) =

--2

0

(1--e-2ZS)}Fl(ds)"

{s--l(1--e-ZS)}F~(ds)

{s--l(1--e-~)}ZF~(ds)0

Z

/~, = g(/3p)

o*cg)o(g) /32

0

+ ph(~3p)W"W

Z

(9)

Theorem. As p -+ 0, while 5 = ~/a and o*(U)o(U)/3 2 remain fixed, the distribution, p , , of the stochastic process: (p)~/2

{

X

(~)

t } --~ ~ ,

0
converges weakly to the Wiener measure with 0 drift vector and dispersion matrix 2D (U), where:

2~(U) = g(3p)o* (U)o(U)/132

(10)

The theorem above may be understood physically as follows. Since the parameter p-i denotes the frequency of collisions of a solute molecule with the solid phase, the theorem states that, when observed at time points far apart eompared to p, X ( ' ) is a p p r o x i m a ~ l y a Brownian m o t i o n with drift vector ~(U)t and the covariance matrix 2D(U)t. The multiplication by pvz (in the theorem) of [X(t/p) -- (t/p)~] merely prevents it from "blowing up", i.e. brings it back to scale as.p-+ 0. The statement in the theorem requiring that as p ~ 0 while 5 and o (U)o(U)/132 remain fixed, amounts physically to an appropriate choice of space and time scales with respect to the parameter p. This aspect is discussed in Bhattacharya and Gupta (1979, p. 487). Note t h a t in the expression for ~ in eq. 9, if one drops the second term ph (3p)W" W, the expression for 2/9(U) is obtained. In order that the theorem be of physical significance, this neglected term in eq. 9 should be very small compared to g(5)o* ( U)o(U)/13 2.

364 T o d e m o n s t r a t e the smallness o f p h ( 6 ) W " W in c o m p a r i s o n with a*(b~)o (U)/fl 2 , a typical range f o r the values o f p m u s t be calculated. Under the cond i t i o n t h a t U = O, a 2 (U =: O)/2fl 2 r e d u c e s t o t h e Einstein's expression for the m o l e c u l a r d i f f u s i o n tensor, dol, f o r a solute in a bulk liquid, where I is a 3 × 3 i d e n t i t y m a t r i x (the n e x t section gives the details o f the a r g u m e n t lead ing to this result). F r o m e x p e r i m e n t a l observations, the ratio of the coef! ficient o f m o l e c u l a r d i f f u s i o n in a p o r o u s m e d i u m , do, to t h a t in a bulk liquid, do, is in the range o f 0 . 4 - - 0 . 8 (Fried and C o m b a r n o u s , 1971). But, f r o m eq. 16 d'o/d o -- g(flp) < 1, where g(flp) is given by eq. 9. As an example, it m a y be assumed t h a t F 1 (s) has an e x p o n e n t i a l density and the expression f o r g([Jp) m a y be e m p l o y e d t o calculate tip t o be in the range o f 2.0--8.0. f o r 0.4 ~ do/d o ~ 0.8. Since typical values o f fi f o r h y d r a t e d ions in electro. lyre solutions have an o r d e r o f m a g n i t u d e in the vicinity o f 10 ~3 s -~' (Robin. son and Stokes, 1965, pp. 43 and 44, and a p p e n d i x 6.1) and 10 '~ ~101° ,'~ f o r t y p i c a l colloidal m o l e c u l e s (see, e.g., Chandrasekhar, 1 9 4 3 , p. 49), it follows t h a t the o r d e r o f m a g n i t u d e o f p is b e t w e e n 10 -9 and 10 " ~~ s. Fo~,' U ve 0, o n e w o u l d n o t e x p e c t p t o be larger t h a n its value at U = 0. It can be easily c h e c k e d f r o m eq. 9 t h a t ph([Jp){W"W) ~ 2 p W " W ~ 2 " 1 0 ~ c m : ',, w h e n [ U] = 0.1 cm/s. F o r this value o f t U 1, the e x p e r i m e n t a l l y observed values o f t h e dmgonal ' e l e m e n t s o f g(fip)a * ( ~U)a( U)/2fi ~ 2 , Le. " the coefflcmnt~ " " o f longitudinal and transverse dispersion, are respectively 10 2 d to and 8d~ (see the n e x t section f o r details). Since do is t y p i c a l l y a b o u t 10 -(' cm 2/s m its o r d e r o f m a g n i t u d e , it is clear t h a t the first t e r m on the right-hand side of eq. 9 is significantly larger t h a n t h e s e c o n d term. Since {X(t); t > 0} has been s h o w n to be a Wiener process, its transitio~ p r o b a b i l i t y d e n s i t y , p(x,t; Xo,t0) satisfies a parabolic differential e q u a t i o n {see, e.g., G i k h m a n and S k o r o h o d , 1 9 6 9 ) , where p ( x , t ; x o , t o ) d v is the probability t h a t a solute m o l e c u l e lies in the e l e m e n t a l m a c r o s c o p i c w~lume du a r o u n d the v e c t o r x at time t, given t h a t at t o it was at x 0 . Since the soluU~ c o n c e n t r a t i o n is assumed to be low, the v e l o c i t y processes {V(t); t ~ 0} cor~ r e s p o n d i n g to d i f f e r e n t solute molecules m a y be c o n s i d e r e d t o be statisti. cally i n d e p e n d e n t . H e n c e , t h e p o s i t i o n processes {X(t); t > 0 j c o r r e s p o n d i n g t o d i f f e r e n t m o l e c u l e s are also i n d e p e n d e n t . L e t the m a c r o s c o p i c solute conc e n t r a t i o n , C(x,t) be d e f i n e d as p r o p o r t i o n a l t o t h e n u m b e r o f solute molecules at t h e p o i n t x at t i m e t p e r unit v o l u m e o f the p o r o u s m e d i u m . T h t ~ law o f large n u m b e r s t h e n can be i n v o k e d t o s h o w t h a t the solute c o n c e n . t r a t i o n , C{x,t) is given b y the e q u a t i o n :

C(x,t) = f P ( x , t ; X o , t o ) C o ( x o ) d 3 x o

4

R 3

where C 0 ( x 0 ) d e n o t e s the initial solute c o n c e n t r a t i o n . T h e details o f this arg u m e n t can be f o u n d in B h a t t a c h a r y a et al. (1976). T h e r e it is p r o v e d t h a t C(x,t), given b y eq. 10, satisfies the same parabolic differential e q u a t i o n as the transition p r o b a b i l i t y density. This differential e q u a t i o n m a y be e x pressed:

365 3

(~C

~2

~ Oij

3 {C(x,t)}-- ~

{o~iC(x,t)}

(12)

Eq. 12 is a general transport equation governing the solute concentration in a homogeneous isotropie saturated porous medium under a steady convective velocity U. The explicit expressions for g(U) a n d / 9 ( U ) in terms of physical parameters were given in eqs. 8, 9 and 10, respectively. It may be noted that, in eq. 12, the mean velocity of the solute, g(U), is different from the mean velocity of the liquid, U. We shall now examine in detail the derived expressions for these two transport coefficients, ~(U) and D(U) and discuss their physical implications.

STRUCTURE OF THE TRANSPORT COEFFICIENTS The structure of the dispersion matrix, D(U), is discussed first under the condition that U = 0. This is followed by a discussion on its structure for a non-zero U. The nature of the convective velocity vector, ~(U), is discussed next in this section. Under the condition U = 0, eq. 4 reduces to the classical Langevin equation (Chandrasekhar, 1943; Nelson, 1967). The "friction coefficient" m~ appearing in eq. 4 then may be assumed to be the same as the one appearing in Einstein's well-known expression for the coefficient of molecular diffusion, do, in aqueous solutions:

do

=

I~ B

T/m~

(13)

where m denotes the effective mass of a solute molecule, ~B is Boltzmann's constant, and T is the absolute temperature. In a bulk solution, o(U = 0), which appears in the Langevin equation, is given by OoI, where 02 ~ 0 and i is a 3 × 3 identity matrix. Then, from the derivation given by Nelson (1967, pp. 54 and 55), one may write 0~/2~ 2 = d 0. Hence, o5 = (2~B T~/m). For a given solute, m~ can be measured directly by an experiment on molecular diffusion in a bulk solution (Robinson and Stokes, 1965). If the solute molecules comprise a large number of atoms (e.g., colloids and complex organic solutes), then it seems reasonable to view these molecules approximately as spheres of radius r. For such situations, a very good approximation to m~ is 6nr~ where v is the viscosity of the solution (This approximation, formally coming through Stokes' law, was also invoked by Einstein.) It has been established from molecular considerations in recent years (Mazur and Oppenheim, 1970). Experimental observations show t h a t the approximation, m~ -~ 6nrv holds reasonably well even for small ions, where r is now the radius of the h y d r a t e d ion (Wolynes, 1978). The typically-estimated values of i~-' in bulk solutions for large molecules and monatomic ions, would be in the range of 10 -9 - - 1 0 -13 S. Neglecting g~-i i n eq. 8, one has ~ ~ 0 and the molecular diffusion matrix, given by eq. 10, D(U = 0) =/90 reduces to:

366

[

1

i ¸t2(1 .... e -~"s) - ~ (11 - e - 2 ~ ' S ) t d F l ( s ) 1I

= dot

(14)

One can now prove that d o ~ d0, i.e. that under zero convection the coefficient of molecular diffusion in a porous medium is less than the coefficient of molecular diffusion in the bulk solution. This inequality is supported by experiments (Fried and Combarnous, 1971, p. 193). In order to prove do ": do, it is enough to show that the factor: j { 2 ( 1 - - e -''~)

~-(1.-e-2'"~)tdF~(s)~0

t151

0

This is straightforward to do, since the integrand in eq. 15 is strictly positive for all s ~ 0; for, it is zero at. s :--- 0 and has a positive derivative for all s ~ 0. A numerical estimation of d o requires that F 1(s) be specified. As before, let us assume that dF 1(s) -= exp(-- s)ds. Then eq. 14 provides an explicit ex! t pression for d o as a function of d0, ~3 and /2. If we take the ratio d o / d o -: 0.67, then for this example, eq. 14 gives 8 = fi/2 2 4. Strictly speaking,/2 -j the mean collision frequency of solute with the solid phase, should be estb mated independently of the measurements of d 0. Independent experimental estimations of/2- ~,/3, d o and d o can provide additional support for the above theoretical development other than what is now available. However, the experimental determination of /2-~ and its functional dependence on other measurable physical parameters (e.g., porosity} constitute areas of further research. Now let I U I ~ 0. Choose the coordinate axes so that U = (U 1,0,0), Uj i: 0. In view of the isotropy of the medium, the dynamics of solute transport, contained in eq. 12, remains invariant under reflection of the coordinate axes. The following discussion will be based only on this assumption of isot r o p y and is thus independent of the molecular m o d e l d e v e l o p m e n t s given earlier. Now one can prove that the dispersion matrix D(U) is diagonal. The proof is based on the argument that if the x2-axis is reflected (replaced by --x 2), keeping the x l-and x3-axes fixed, the covariance terms involving x 2 change sign. From the definition of isotropy, such a reflection does not alter D(U). Hence, D12 ....... Dl2 :-= 0, D23 = --D23 = 0. Similarly, D13 :-- 0. Since D(U) is symmetric, it follows that it is diagonal. It does n o t follow from these arguments t h a t the three diagonal elements of the dispersion coefficient matrix D(U) are e q u ~ . Indeed, it has been observed experimentally that the diagonal elements of D(U) increase as functions of U and that D 11 is "always larger than D22 and D33 (Fried and Combarnous, 1971). With this in view, we can expand the Dii as functions of U in a Taylor series around U = 0. The first term in the expansion is clearly d o . The expansion may be written in the form :

367

Dii = do +ailU1 + a n U ]

+...,

i = 1,2,3

(16)

It now follows from the invariance of Dii, i = 1, 2, 3, under the reflection of the xl-axis , that Dii must be even functions of U~. Thus eq. 16 reduces to:

Dii~--do +ai2U~,

i = 1,2,3

(17)

According to eq. 17, the dispersion coefficients are approximately parabolic as functions of U 1 for small U = (U1,0,0). Finally, it may be noted that isotropicity in x2--x 3-plane implies that D22 = D 33- A physical specification of the coefficients a12 , a22, etc., appearing in eq. 17 requires further research on both the theoretical and the experimental aspects of solute transport in porous media. The above discussion on the structure of/9~U), for I U I =# 0, was independent of the molecular model expression for D(U), developed in this paper. Now, on the basis of the expression for D(U), given by eq. 10, we will show t h a t the molecular parameters o ( ' ) and/3(') in the Langevin equation must depend on the mean convective fluid velocity U. In order to see this, we first observe that, as per the proof given after eq. 15, the function g([J#) ~ 1, irrespective of the values which ~p takes. However, experimental results show that for moderate values of [U I, the coefficients of longitudinal and transverse dispersion, i.e. D l l and D~2, are much larger than d~; their c o m m o n value for I U I = 0. As an example, if we take the Peclet number in Fried and Combarnous (1971, p.228) to be [ U ['103 , then for I U[ = 0.1cm/s, the coefficient of longitudinal dispersion, D l l , is 102d0 and the coefficient of transverse dispersion is 8do. Hence it follows from eq. 10 that o and ~ in eq. 10 must vary significantly with the mean convective velocity of the fluid. Therefore, the dependence o f / ) ( U ) on U primarily comes through the dependence of a and ~ on U, because g(fip) is always less than one. This suggests that the dispersion effect, manifested by the dependence of 9 ( ' ) on U, primarily comes through the solute--liquid collisions rather than through the collisions of the solute molecules with the solid phase. However, the dynamical reasons underlying the appearance of the dependence of o and ~ on U are not understood at present. An investigation of this aspect constitutes a very important area of basic research in the realm of solute dispersion in porous media, because it lies at the root of what is c o m m o n l y referred to as dispersion in porous media. Now, in order to understand the structure of the mean convective velocity of solute, ~(U), consider the effect of 5 = ~p on it in eq. 8. It can be checked simply that (see eq. 8):

limf(5) = lim [i -1~- ;(l--e-~)dF1(s)l--+l 5t~ 5t~ o

(18)

Therefore, when ~p is large, it follows from eqs. 18 and 8 that ~ ~ U, i.e. the solute convective velocity is approximately equal to the convective velocity

368 of the liquid. On the basis of experimental observations, it has been often stated in the literature, that ~ --~ U (see, e.g., Nielson et al., 1972, p.123), especially for unconsolidated media, e.g., glass beads, sands, etc. The first conclusion that can be drawn from the observed equality of a and {f, at least for some media, is t h a t 5 = fig must increase as a function of the mean liquid velocity U. Let us demonstrate this through an example. As before, assume t h a t dF(s) = exp(-- s)ds. Then the computations following eq. 15 show that for U = 0, and for a typical value of the ratio do/do = 0.67, t # = 4. For this value of tip, f(flP), given by eq. 8, is 0.8. This result, in conjunction with the observed equality of ~ and U for I UI =~ 0, implies that tip depends on U and is larger for non-zero U than when U = 0. At this stage it is not (:lear as to what extent t # will vary with U for any arbitrary homogeneous porous medium. For example, will ~ and U be the same even for unconsolidated media? A full understanding of the behavior of the mean convective solute velocity awaits further research on both the theoretical as well as the experimental aspects of solute movement in porous media.

SUMMARY AND CONCLUSIONS The principal results in this paper may be summarized as follows: (1) A molecular model of solute transport has been developed from three stochastic postulates concerning the dynamics of a solute molecule traversing a porous medium. (2) A rigorous mathematical consequence of these postulates is a parabolic differential equation for solute transport in water-saturated homogeneous isotropic porous media, as well as expressions for the convective velocity vector of solute and for the molecular diffusivity in terms of physicallymeasurable parameters. (3) The convective velocity vector of the solute in a porous medium is different from (less than) the convective velocity of the liquid in a porous medium. Theoretically speaking, the convective solute velocity can become approximately equal to the mean liquid velocity. To what extent this difference is significant in various porous media is not well understood at present. (4) The coefficient of molecular diffusion in a saturated porous medium, under zero convection, is shown rigorously to be always less than the coefficient of molecular diffusion in an aqueous solution. Experiments confirm this deduction. (5) S y m m e t r y arguments are developed to show that, if the c o n v ~ t i v e velocity vector U = (U I ,0,0) is oriented along the principal axes, then the dispersion matrix is always diagonal for an isotropic medium. Because the s o l u t e dynamics, given by a parabolic equation, can be assumed to be invariant under reflection of the coordinate axes if the medium is isotropic. An additional consequence of this invariance is that for small U 1 , the behaviour of the dispersion coefficient as a function of U~ is approximately parabolic.

369 (6) It is suggested t h a t w h a t is c o m m o n l y r e f e r r e d t o as dispersion in a p o r o u s m e d i u m c o m e s primarily t h r o u g h the collisions o f the solute with the liquid m o l e c u l e s r a t h e r t h a n t h r o u g h the solute--solid phase collisions. However, the d y n a m i c a l basis for the a p p e a r a n c e of this b e h a v i o u r is n o t u n d e r s t o o d at present. It is n o t clear as to w h a t role the p o r o u s m e d i u m plays in p r o d u c i n g this d e p e n d e n c e . F o r e x a m p l e , is this d e p e n d e n c e present even in liquids flowing slowly t h r o u g h capillaries? Does it p o i n t to a basis on which the t h e o r y d e v e l o p e d b y T a y l o r ( 1 9 5 3 ) needs r e e x a m i n a t i o n ? An answer to this q u e s t i o n and m a n y o t h e r i m p o r t a n t questions, e.g., h o w d o the results o f this w o r k e x t e n d t o immiscible solutes, and h o w can one i n c o r p o r a t e the behavior of precipitating and adsorbing solutes into the m o l e c u l a r d y n a m i c s , await f u r t h e r research. In a r e c e n t review article ( S p o s i t o et al., 1979), it was m e n t i o n e d t h a t the c o n c e p t u a l t h i n k i n g b e h i n d the h y d r o d y n a m i c a p p r o a c h e s to the f o u n d a t i o n s o f solute t r a n s p o r t t h e o r y m a y n o t be physically and m a t h e m a t i c a l l y justifiable. This p a p e r gives q u a n t i t a t i v e insight along these lines. S u p p o s e U = 0. In this case, the parabolic diffusion e q u a t i o n in bulk s o l u t i o n is valid o n l y at a time scale which is very large c o m p a r e d to ~-1 (see, e.g., Chandrasekhar, 1 9 4 3 ; Nelson, 1967). Thus, in o r d e r t o invoke a diffusion a p p r o x i m a t i o n at the p o r e level (a f u n d a m e n t a l a s s u m p t i o n b e h i n d the h y d r o d y n a m i c approaches), it is necessary f o r a solute m o l e c u l e to remain in the liquid phase for times which are very large in c o m p a r i s o n with ~-~. H o w e v e r , o u r c o m p u tations f o l l o w i n g eq. 10 show t h a t ~ 1 and p are of the same o r d e r o f magnitude. Thus the p r e s e n c e o f the solid phase prevents a solute m o l e c u l e f r o m traveling u n i m p e d e d in the liquid phase for times which are very large comp a r e d t o ~-1. This shows t h a t the validity o f t h e diffusion a p p r o x i m a t i o n at the p o r e level is on d u b i o u s grounds.

ACKNOWLEDGEMENTS T h e research r e p o r t e d in this p a p e r was s u p p o r t e d b y Grants CME 8 0 0 4 4 9 8 , CME 8 0 0 4 4 9 9 and CME 7 9 2 0 7 7 8 f r o m the National Science F o u n d a t i o n .

REFERENCES Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York, N.Y. Beran, M.J., 1968. Statistical Continuum Theories. Wiley, New York, N.Y. Bhattacharya, R.N., 1978. Lecture notes on stochastic processes. Math. Dep., Univ. of Arizona, Tucson, Ariz. Bhattacharya, R.N. and Gupta, V.K., 1979. On a statistical theory of solute transports in porous media. Soc. Ind. Appl. Math., J. App. Math., 37 : 485--498. Bhattacharya, R.N., Gupta, V.K. and Sposito, G., 1976. On the stochastic foundations of the theory of water flow through unsaturated soil. Water Resour. Res., 12 : 503--512. Chandrasekhar, S., 1943. Stochastic problems in physics and astronomy. Rev. Mod. Phys., 15 : 1--89.

370 Chu, S.Y. and Sposito, G., 1980. A derivation of the macroscopic solute transport equation for homogeneous, saturated, porous media. Water Resour. Res., 16: 542--546. Fried, J.J. and Combarnous, M.A., 1971. Dispersion in porous media. Adv. Hydrosci., 7 : 169--282. Gikhman, I.I., and Skorohod, A N . , 1969. Introduction to the Theory of Random Pro ~ cesses. W.B. Saunders, Philadelphia, Pa. Mazur, P. and Oppenheim, t., 1970. Molecular theory of Brownian motion. Physica, 50: 259--288. Nelson, E., 1967. Dynamical Theories of Brownian Motion. Princeton University Press, Princeton, N.J. Nielson, D.R., Jackson, R.D., Cary, J.W. and Evans, D.D., 1972. Soil Water. Am. Soc. Agron., Madison, Wis. Present, R.D., 1958. Kinetic Theory of Gases. McGraw-Hill, New York, N .Y. Rice, S.A. and Gray, P., 1965. The Statistical Mechanics of Simple Liquids. Wiley~ New York, N.Y. Robinson, R.A. and Stokes, R.H.~ 1965. Electrolyte Solutions. Butterworths, London. Scheidegger, A.E., 1974. Physics of Flow Through Porous Media. University of Toronto Press, Toronto, Ont. Sposito, G., Gupta, V.K. and Bhattacharya, R.N., 1979. Foundational theories of solute transport in porous media: a critical review. Adv. Water Resour., 2: 59--68. Taylor, G.I., 1953. Dispersion of soluble matter in a solvent flowing slowly through a tube. Proc. R. Soc., London, Ser. A, 219: 186--203. Todorovic, P., 1970. A stochastic model of longitudinal diffusion in porous media. Water Resour. Res., 6: 211--222. Todorovic, P., 1975. A stochastic model of dispersion in a porous medium. Water Resour. Res., 11: 348--354. Uhlenbeck, G.E. and Ornstein, L.S., 1930. On the theory of the Brownian mo~ion. Phys. Rev., 36: 823--841. Wolynes, P.G., 1978. Molecular theory of solvated ion dynamics. J. Chem. Phys., 68: 473--483.