- Email: [email protected]
Multicomponent solute transport with sorption and soluble complexation
Multjcomponent solute transport with sorption and soluble complexation DAVID J KIRKNER, THOMAS L. THEIS and AARON A JENNINGS Department of Ctvtl Engmeermg, Umverstty of Notre Dame, Notre Dame, IndLana 46556, USA
A numerical algonthm for predicting the migration of multiple subsurface pollutants has been developed accounting for dispersion, convection, soluble complexatlon and sohd phase accumulations (sorption) The basis of the model is a fimte element solution of the mass transport equation The essence of the algonthrn is the treatment of the sorptlon terms as implicit functions of the total soluble concentrations facilitating a general and modularized treatment of solution phase and solid phase chemistry Examples are presented dlustrating the effects of soluble complexatlon and competitive sorptlon on the transport of multIcomponent solutions
INTRODUCTION The movement of groundwater contaminants has received increased attention In recent times Both organic and inorganic pollutants from such sources as solid wastes and sludges disposed of on land, wastewaters renovated through groundwater recharge, and polluted surface waters which recharge groundwater aquifers are examples where the need for an understanding of the pnnclples governing mass transport exists Frequently complex physical, chemical, and biological interactions are Involved and must be included in the analysis of the problem. The purpose of this paper is to present a numerical algorithm to solve the equations governing the migration of subsurface pollutants accounting for dispersion, convection, soluble complexatlon and solid phase accumulations (sorptlon) Several researchers have addressed the problem of the transport of chemically reacting substances employing different methods and with varying degrees of generality For multlcomponent transport the problem is, In general, governed by a system of nonhnear partial differential equations This system of equations can be formulated explicitly If the relationships between the solid and the soluble form of the species are such that the solid phase terms can be eliminated from the transport equations Rubin and James ~ formulated the governing equations for multlcomponent one-dimensional, transport where ion exchange ~s the sorptwe mechanism The resulting equations were solved by a Galerkln finite element method (FEM) using a Crank-Nlcolson predlctor-corrector time marching scheme Valocctu et al 2 extended this technique to two-dtmenslonal domains Jennings Klrkner, and Thels 3 formulated the govermng equations for multlcomponent, one-dimensional transport Including soluble complexatIon and nonlinear adsorption These equations were solved using the FEM with a backward difference time marching scheme and a Newton-Raphson iteration within each tinle step Lln 4 formulated the problem for one-dimensional transport in a two-layered medium Including nonhnear adsorption with a Freundllch Isotherm He reduced the PDEs to ODEs by orthogonal collocation and then rotegrated usIng a fourth-order Runge-Kutta scheme Accepted April 1984 Discussion closes November 1984
120
Adv Water Resources. 1984, Volume 7, September
In all cases, above, substituting the chemical equations into the transport equations causes some or all terms to be nonhnear An alternative approach IS to omit the substitution and to treat the problem as mass transport wath a nonlinear source term This approach was taken by Trotta and Del Giudice s and Winters et al 6 in modeling coupled heat and mass transfer If the numerical scheme is such that the source term must only be evaluated but need not be known explicitly, then complicated chemical systems can be handled. In this paper the latter techmque is adopted The essence of this algorithm is the treatment of sorptive terms as tmphclt functions of the total soluble concentrations It will be shown that th~s allows a very general and modularized treatment of solution phase and sohd phase chemistry
NOTATION Bold face mlnlScules u, ~, etc represent vectors The number of components of each vector wdl be given where first used or else should be clear from the context Bold face maglscules are matnces K, M, etc The dot product of two matrices of the same order is defined A-B = Y,,IA,IBu A superposed dot indicates a spatial time derivative, u = (O/Ot)u(x, t) The gradient of a vector is the matrix Vu with components (Vu)u = Ou~/Ox/ The divergence of a vector u is defined dlv ( u ) = Z, Ou,/Ox, Lightface symbols are scalars u,, O, t, Atr etc
STATEMENT OF THE PROBLEM To help set the framework for this discussion, consider the simple three-component system with 1 2 ion exchange and soluble complexatlon, governed by the following reactions
2Cl + 62 ~ 2 ~ + c2, Ks cl +c3 ~ c l 3 , K 1
(1)
C2 + C3 ~" C23, K2
Cl and c2 are adsorbing metals and c3 is a nonadsorbmg hgand which undergoes solution phase reactions with the metals to form the soluble complexes, c13 and c23 6~ and 62 are the adsorbed forms of cl and c2 respectively For a saturated batch system a mass balance per unit volume for
0309-1708/84/030120 $2 00
© 1984 CML Pubhcanons
Multtcomponent solute transport wtth sorptton and soluble complexatton D J Ktrkner, T L Thels, A A Jennmgs each component yields after substituting the mass action equations from ( 1 )
where 0BI W 0B2 = 0B is the complete boundary and n is the outward unit normal and u(x 0)=U,o(X)
/ e l / + K I [ C 1 ] [c3] -k--~{Cl} -- [CIT]
[e~l +/,'~ [o~1 [c31 + -p0[e:IK_~( [e,] IIe,]i~ = [e~rl
(2)
[e3] +Ki[Cl] [e3] +K2 [c2] [e3] = [e3T] where [ ] = mass/volume of solution, { } = mass/mass of solids, and [c,T ] is the specified total concentration of component t Here p is the bulk density and 0 is the porosity of the domain The remmmng required equatmn is the mathematical expression of the definition of the exchange capacity, { dT}, 1 e
(¢1} + {e:) = {er}
(3)
or using the mass action equation from (1),
{if'l} + ~
\~'~!
= {CT)
(4)
Equations (2) and (4) are four equations for the four unk.nowns [cl], [c2], [ca], and [61] (they could alternatively be expressed with [c2] as a pnmary unknown) The appropriate mass achon equations then allow the soluble complex concentrations, [c13] and [c23] and the remmning sohd, {dz} to be determined For a flowing system equation (2) is simply a statement of the fact that the total concentration of each component is the sum of the species solution and sohd phase concentrations The conservation of mass equations for a flowing system must additionally be written This is facdltated by the introduction of the total soluble concentrations as variables, 1 e let (from hereon the brackets denoting concentratmn will be omitted) U 1= C1
+Kxcxe3
U2 = C2 ~/k'2c2e 3
(5)
u3= e3 + KIelC3 + K2c2e3 Now equation (2) can be rewritten I'll
+_o
(1 = ClT
.{_._PPC2 ( ( 1 ) 2 "c-~1" = c2v
0 KS
(9)
where B IS a hounded regular region Equations (5)-(9) constltute a complete statement of the problem [inphclt in the above is the local equilibrium assumption, that is the reactions, (1), are assumed fast enough compared to the transport process that they may be assumed instantaneously at equilibrium. The O~/Ot term an (7) couples the differential equations and makes them nonlinear in general This virtually rules out analytical solutions There are, however, inany approaches possible for a numencal solution Jenmngs et al 3 subsntuted equanon (5) Into (7) and also expressed (: in terms of el, cz, c3 (they considered Langmulr type adsorpnon) thus converting the governing differential equations with u as dependent variable to a differential equation with species concentration e = (c~, e2, e3) T as dependent variable This approach was also used by Valocchl et al 2 for an ion exchange system similar to that considered here only without soluble complexation The disadvantage of this technique IS that the program to solve the differential equations becomes very problem dependent This makes changing the number of components or the form of the solid phase partinonlng difficult Additionally, as the complexatlon increases (1 e as K1 and K2 become large forcing a greater percentage of the mass to exist as a nonadsorbmg complex), the solution of the differential equations becomes Increasingly difficult 3 The techmque presented here retains u as the primary dependent variable In (7) and treats ~ as an lmpllot funcnon o f u This method, as will be shown, keeps the chemical interaction calculations separate from the transport calculations and in fact allows the solution phase chemistry to be treated separately from the solid phase chemistry The method of dlscretazatlon Is the Galerkln fimte element method In space and a fimte difference method m time The resulting nonhnear algebraic equanons are solved by a Newton-Raphson procedure Other methods for discretizatlon could be employed What is emphasized here is the treatment of the chemical mteracnons due to the sorpnon term xn (7) SOLUTION ALGORITHM
0
u:
xEB
(6)
Taking the dot product of both sides of (7) with an arbitrary vector w and integrating over the domain, using the dwergence theorem and (8) yields the weak form of (7)
U 3 = C3T
The conservanon of mass equanons for saturated, lsochonc, groundwater flow with sorpnon are 7 dw (D Vu) = (Vu)v + u + (p/O) e
(7)
where u IS the vector of total soluble concentrations, 5 are the sohd concentrations, D IS the dlspersmn tensor, and v is the interstitial velomty field (dlvv = 0) For the sample problem introduced above
f w'qdA
B
i~B:
(10) where w E H ° = { w l f ( w ' w + Vw'Vw) d V < ~', w(x) = 0, x E 0BI}
B Construct finite element Interpolations for approximate solutions to (10) as follows Let
I1 = (Ul, U2, U3) T e =
f {(DVu)'Vw+w'[(Vu)v+u+(p/O)~]}dV=
(6~, d2, O)r
N'M
Statable boundary and lmtlal condltmns are reqmred with (7), for example u(x, t) = fi(x, t), x E 0BI D(Ve) n = q, x E 0B2
(8)
U ~-u h =
~
~laj(t)
1--=1 N,M
(11)
(p/O)e ~- (o/O)e~ = ~ ~lbl(t) 1=1
Adv WaterResourees, 1984, Volume 7, September
121
Multtcomponent solute transport with sorptton and soluble complexatton D J Ktrkner, T L Thets, A A Jennmgs where N is the number of nodal points and M iS the number of chemical components (in the sample problem above. M = 3) In equation (11)
¢~ : (~i. ~ ,
, ~m) T
(12)
are finite elemei~ basis functions with the following property ~]k(X~)=81t,
l=k+(t--1)M
(13)
where x~ means the tth nodal point and 6tt is the Kronecker delta Equation (13) lmphes
degrees of freedom in the standard manner Finally (17) can be wntten m vector notation as follows Ka + Ma + Mb = f
The global matrices m (19) may be assembled from element contributions m the standard manner s Equation (19) is called the semldiscrete Galerkm approximation Discretlzatton m tune further reduces the set of ordinary differential equations, (19), to algebrmc equations A vanable two-step time marching scheme is employed here s Let 1 a=--(an+l--an), At
Uk(X,) = at
(p/O)Ck(X,)=bt,
(14)
l=k +O--1)M
Thus the coefficients in the expansions (11) are nodal values of soluble concentrations and adsorbed solid concentrations This means also that only one term m each vector ~/is nonzero This is expressed mathematically as $jk(x) =
{
:~ 0, = 0.
k = l(mod M) k :/:l(modM)
(15)
N'M I [(DV~j)'V~, + ~, "(V~/v) ] dVa] • B
+
(c#,.¢])dVb]
(¢,-¢j)dVa t + ./:1
LB
* B
=
(¢,-q)dA aB:
(16) or
N.M (K,/a I +Mqa I +M,tb/) =f~,
t ~ [ t n , tn+l]
1
(20)
b=--(bn+|--bn), At
t E [ t n , tn+,]
and
a(tn÷cO = aan+l + (1 - - a ) a n
l = 1,N'M
(17)
/=i
where K,j = i [V$,'(DV~bj) + ¢,'((V¢t ) v)] dV
with 0 < a < 1 In the above
tn+ ~ = (n + a) At Evaluating (20) at tn+ ~ yields 1 l(an+' + ~tt Mbn+| = p
l(=aK
1
+--M At
1 p = f + ( I ( - - K) an + Att Mbn Various standard time marching schemes result by a suitable choice of o~ For instance a = 1 is the backward difference scheme and a -- ½ is the Crank-Nlcolson method Equation (23) is a set of nonlinear algebraic equations The simplest lteratwe scheme for solving these is to assume a guess for bn+ 1, say 0, and solve for an÷l Then solve the chemical problem at each node (equations (3), (4) and (5)) to obtain an unproved value for bn+ 1 This procedure IS then repeated until convergence Note that because of equation (15) I( and M have no off diagonal temls that represent coupling between components The only couphng that occurs is between degrees of freedom representing the same component at different nodes Thus (23)can actually be solved as M one component problems The principal &sadvantage of this approach is that It does not always converge In fact for many cases of practical importance this simple method of successwe approximations is madequate An alternative procedure is to solve (23) by a NewtonRaphson lteratwe technique 9 This can be wntten
where
],= f
(18) g = p -- I~an+l
(¢,-q) dA
~)B2 The equatJons (1 7) can be reduced to account for specified
122
Adv Water Resources, 1984, Volume 7, September
(23)
where
m+| rr/ an.l = an+| + (jm)-lgm
B
(22)
an = a(tn)
B
M,] = f (q~,'¢/) dV
(21)
b(tn+~) = abn+l + (1 - - a ) bn
It should be noted that the shape functions for the adsorbed solids do not have to be the same as for the solutes In fact since the adsorbed solids are considered implicit functions of the solutes only the inner product of the test functions and the vector of adsorbed solids has to be evaluated, the actual values of the adsorbed solids are never required There are two reasons for adopting the interpolation scheme of equation (11), however First the amounts of the adsorbed solids are of intrinsic interest Second a substantial amount of computer time can be spent evaluating the inner product mentioned above, since it will be shown to involve chemical speciatlon When the solids are interpolated as m (1 I) the number of chemical calculations is kept moderate whde also reducing the inner product calculations since the M matrix (see below) only needs to be computed once Now using (1 1) m (10) and taking the test functions w to be ~ (t = 1, N . M ) yields
/=1
(19)
1 ~tt Mbn+|
1 J : r, + - -
At
MJ
J = (ab/Oa) [tn+,~
(24)
Multwomponent solute transport with sorptlon and soluble eomplexatton D J Ktrkner T L Thets, A A Jennmgs and the superscript m ln&cates the iteration number The principal difficulty w~th this method hes m the calculation of the Jacoblan matrix, J, since b is given as a function of a only implicitly Note that by virtue of equation (15), J has a block dmgonal form (see Fig 1) where each submatnx is itself a Jacobmn matrix of the form Jt = (Oe/OU) at (xt, tn.~),
l = 1,N
(25)
and each J~ is of dimension M x M Since e is an implicit funcnon of u, one approach for evaluating Jt is to calculate the derivatives numerically lo Oat
at(Ut,
, U1 + 6,
,Um)
0u /
-
-
Ct(U 1,
, U 1,
, UM )
Thus e = goh-l(u)
(30)
(ae/0e) Oe/Ou)
(3 l)
and then from the chain rule
~e/0u
=
Since equanon (28) is often solved numerically m chemical equ~hbrlum programs by Newton's method, the Jacoblan matrix 0h/0c is available lz From equanon (28)
ac/0u =
(0h/be)-'
(32)
For many exchange expressions (such as the mulncomponent Langmulr isotherm), the matrix 0e/0e may be evaluated explicitly Thus equation (31) may be written
e
Oe/au = (~e/Oe) (0h/0e) -t
where e is some small parameter To evaluate each term on the right requires (a) Solution phase speclanon For the sample chemistry introduced previously th~s means solving equation (5) for ct. c:, c3 gwen u (b) Sohd accumulation Given c~, c2, cs, use equation (4) and the mass action equation for (1) to calculate 61 and (2 The key aspect of the above algorithm ~s that the chemical computations may be modularized, 1 e kept ennrely separate from transport calculanons In fact solution phase computations are distract from sohd phase computations thus facilitating the study of different sorptlve mechanisms Some degree of generahty may be obtained ff the solunon phase specmnon calculations (part (a) above) are performed by one of the available general-purpose chemical equthbna programs (e g Morel and Morgan H) It must be noted that when a fairly comphcated chemical module is employed the computanonal rime reqmred to form the Jacobmn matrices, Jt may become rather burdensome Use of a modified Newton iteration, I e use jo instead of j m in equation (24). reduces computatmnal nine considerably and has been found to work well An alternative method, however, may be employed to calculate Jt by using the reformation avadable m the solunon phase chenucal module First, to generahze equahon (5) define the vector of species concentranons as
c = (~ l,
= g(h-l{u))
, CM)T
(27)
Then equanon (5) can be written u = h(e)
(28)
where h Is m general a nonhnear operator Equation (4) and the mass action equatnm tar (1)~ can be generahzed to = g(c)
Equatmn (33) may be used in lieu of equation (26) Equation (26) requires M + 1 calls to the chemical module of the program to evaluate J . whereas equatmn (33) only requires one call However, (33) requires the inverse of 0h/0c The reverse of 0h/0c nonetheless is easier to compute than speclatmg M times Also note that equations (33) or (26) are best computed globally and thus the global matrices, K and M need only be computed once Results of a sample problem wtll now be presented to demonstrate the type of problem that can be addressed with th~s algorithm
EXAMPLE
In ref 1 Rubm and James studied several systems where ion-exchange was the sorptlve mechanism One particular example involved the leaching of a sodlc soft governed by a 1 2 binary exchange reaction The transport is one-dimensional with a mixed or tlurd type boundary condmon at the source and a zero gradient at the downstream boundary (to s~mulate a semHnfimte region) This hypothencal problem wdl be reexamined here and the effect of the addmon of a nonadsorbmg hgand (SO~) which forms soluble complexes with the adsorbing metals will be stu&ed It should be emphasized that this example is not offered as a thorough analysis of the Ca/Na/SO4 transport system, but rather to dlustrate the potential of the method presented hereto to study such systems Thus the governing reactions are equation (1) with cl = Na +, ~2 = Ca m and c3 = SO4 The mass balance and boundary conditions for thts problem are 02tl~
OU t
D--=v Oxz
Ox
+u,+(p/O)6,,
t=1,3(¢3=0)
(34)
I = 1,3
(35)
(29)
( --D'~--x Ou, + vu,} = v~,, "lx=O Out
DEGREE OF FREEDOM NODE
1
I I
2
f I
. . . . . . .
(33)
i I
x=r=0'
1=1.3
(36)
N
The m~ttal condlhons are
co--iA = Jl,,,oEJ21-i,.., o., , liON] Ftgure 1 Contrtbutton to the Jacobtan matrix oJ the sorptton terms
ut(x,O)=uto,
1=1,3
(37)
Numerical values for all the physical ~.onstants, boundary condmons and chenncal data are given m Table 1 Rubm and Janles considered three cases for this problem The firs! one (ref 1 case H) had a high calcmm lnfluent and no sulfate (U3o = u3 = 0 0) The second case (ret 1 case L) was identical in all respects to the first case except the calcmm mfluenl concentration was low The
Adv Water Resources, 1984, Volume 7, September
123
M u l t w o m p o n e n t solute transport with sorptton and soluble complexatton D J Ktrkner, T L Thezs, A A Jenmngs Table 1
Datafor sample problems
160
ut
ut0
Component
1 2 3
0 0020 0 0072 0 100
0 0280 0 0020 00
Na Ca SO,
1
I
I
I
140
Boundary and amtml conditions (meq/cm 3) l
I
120 100
u,
80
/ /A / /...'o.,
G, 6 0 40
Physical data Symbol
Value
Umts
0
0 54 1.20 25 36 2 536
g/cm 3 cm2/day cm/day
....
P'~ '--~-
D v Chemical data Symbol
Value
K, K: KS
2 22 79 2 0 078
o
1' 0
I
I
40
80
1
I
I
120 160 D E P T H , cm
7'00
240
Figure 2 Normahzed sodium concentration versus depth Ul - total soluble sodium concentration, u2 - total soluble sulfate concentration
Umts cm3/meq * cm3/meq * cm 3/gt
I0 U z w, t h ~ 3
~
08 '~
* Equfl]brmm constants trom ref 13, corrected to an average constant xomc strength of 0 05 using the Davies equation From ref 1
\
third case was similar to the second case except gypsum was present in the soft It was assumed that the gypsum was present in a sufficient amount so that the calcium-sulfate solubility product was always sausfied This allowed the number o f operational equations to be reduced to two Herein the first case is reexamined, but w~th sulfate added to the mfluent source Sulfate has been added at a concentration that will allow for complexaUon, but not precxpltatlon o f calcium sulfate Thus there are still three operational equaUons (equation (34)) This example serves to demonstrate the application of the algorithm presented m this paper to a problem with soluble complexation and adsorption Linear basis funcuons are employed for the fimte element Interpolation o f solutes and solids, and a backward difference time marching scheme is used Figure 2 shows normahzed concentration profiles for total soluble sodium at 14 6 and 43 8 days The sohd lines are results with no sulfate mfluent, agreeing with Rubln and James' case H The dashed lines are the total soluble sodmm profiles when the sulfate mfluent is 0 1 0 m e q / c m 3 The presence o f the sulfate has a small effect on the sodium profiles Figure 3 shows normahzed concentration profiles for total soluble calcium w~th and w~thout sulfate present and for total soluble sulfate The presence o f the sulfate causes a calcmm pulse to propagate thiough the domain along with the sulfate front This effect can be explained as follows Figure 4 shows several batch isotherms for this chemical system, 1 e adsorbed c~ (or c:) versus c2 for fixed c~ F o r the simulation presented above, c~/clo is always less than 14 and thus the soil has a strong affinity for c2 except at very low concentrations o f c2 (Ca 2+) Given the original pore water concentrations the behavior o f the profiles without sulfate is to be expected Ion exchange only affects transport in the mmal portion o f the domain
Adv Water Resources, 1984, Volume 7, September
=00
u z w l l h "~5 = 0 ' 0 .
~
~
'
-
U3
-
\
U, 04
124
0 ~
20 o
p
~u 3 " :
.~.~"
ok]:~--;~~--0
40
,\''"r'--~
,
80
120
160
DEPTH,
200
240
cm
Figure 3 Normahzed sulfate and calcium concentration versus depth u ~ - total soluble sodium concentration, u2 - total soluble calcium concentration, u3 - total soluble sullate concentration
I0-
)
CI /C)o= 5
i
i
I
0
.
08
02
06
04
Ez
T,
CT 04
06
02
08
0
0
I
02
I
04
i
06 C2/C20
I
08
TT
I0 0
Figure 4 Sorbed calcium versus soluble calcium concentration cl - [Na÷], (1 - sorbed sodium, c2 -[Ca++], c2 sorbed calcium, ET - cation exchange capacity
M u l t t c o m p o n e n t solute transport wtth sorptton and soluble complexanon D J Kwkner, T L Thezs, A A Jennmgs
0 8
principal advantage o f this approach is the modularlzanon o f the chemical calculanons Solunon phase calculanons are also performed separately from solid phase calculations Add]nonal one- and two-dimensional mult]component apphcat~ons with complex chemical interactions are m preparation It should also be mentioned that tMs approach readily extends to the problem o f mass transport with sorpnon governed by chemical kmet]cs 14
8 Ooys 146 Doy s - ' ~ cz I u z
c~., un
~ C l /
UI
04
ACKNOWLEDGEMENT 02
0
I
40
80
I
I
120 160 DEPTH, cm
I
I
200
240
2BO
Ftgure 5 Fracture o f each c o m p o n e n t as free ton versus depth c] - [Na÷], u i - total soluble sodtum concentratton, c2 - [Ca*+], u 2 - total soluble calcium concentratton
V,qth sulfate present the calcmm (uz) in solunon does not all exist as Ca 2÷ (c2) and the sodium ( u 0 m solution does not exist only as Na ÷ (cO because o f c o m p l e x a n o n With sulfate present the calcmm (sodium) ex]sts as the exchanging free ion, Ca 2÷ (Na*) and as the complex CaSO~(aq) (NaSO~) Figure 5 shows a plot o f the free ion concentration o f each c o m p o n e n t as a fraction o f the total soluble concentration versus depth at 14 6 and 43 8 days Notice that not much sodmm is complexed However, m regions o f Mgh sulfate concentranon (compare with F~g 3) a large percentage of the calcmm is complexed Recalling the Isotherms in Fig 4 it is expected that adsorbed sodmm decreases while adsorbed calcmm increases Since this ~s 1 2 exchange, one should expect a greater increase m the total soluble sodmm than the corresponding decrease m total soluble calcmm This is indeed the effect displayed m Figs 2 and 3 In the original smmlanon w~thout sulfate, ion exchange only affects the first 10 cm o f the domam With sulfate added, a slgmficant degree o f soluble complexatlon occurs (especially CaSO4°) which affects ion exchange as far into the domain as the active sulfate front has penetrated CONCLUSIONS In this paper, an algorithm has been presented for solving the mass transport equanons with sorpnon and soluble c o m p l e x a n o n The procedure employs a Galerkm finite element spanal d~scretlzanon and a finite difference time marching scheme S o r p n o n terms are treated as lmphclt functions o f the total soluble concentranon variables The
The work presented m this paper was supported by grant DE-AC02-79EV 10253, from the Ecological Research Dwlslon, Office o f Energy Research, US Department of Energy to the University o f Notre Dame
REFERENCES 1 Rubln, J and James, R V Dlsperslon-attected transport of reacting solutes in saturated porous medm Galerkm method apphed to equxhbrmm controlled exchange m umdlrechonal steady water flow, WaterResour Res 1973, 9 (5), 1332 2 Valocchl A J , Street, R L and Roberts, P V Transport ot ion-exchanging solutes m groundwater chromatographic theory and held simulation, Water Resour Res 1981, 17 (5), 1517 3 Jennmgs, A A , Ktrkner, D J and Thels, T L Multi-component equfllbrmm chemistry m groundwater quality models, Water Resour Res 1982, 18 (4), 1089 4 Lln, S H Nonlinear adsorption m layered porousmed]a flow, J Hydrauhcs Dlv, ASCE 1977, 103 (HY9), Proc Paper 13192, 951 5 Trotta, A and Del Gmdlce, S A finite element solution of coupled heat and mass transfer with chemical reaction, Chem Eng Sct 1980, 35, 709 6 Winters, K H , Rae, J , Jackson, C P and Chffe, K A The finite element method for laminar flow with chemical react]on, Int J Num. Meth k.ngng 1981, 17, 239 7 Bear, J Dynamws o f Fluuts m Porous Medta, American Elsevier Co , New York, NY, 1972 8 Zlenk]ew]cz, O C The Fmtte Element Method, McGraw-Hall, 1977 9 Ortega, J M and Rhemboldt, W C ]terattve Solutton of Non hnear Equattons zn Several Vartables, Academic Press, New York, 1970 10 Oden, J T Fmtte Elements of Nonhnear Contmua McGrawHdl, 1972 11 Morel, F and Morgan, J J A numerical method for computing equlhbna m aqueous chemical systems, finvzron S c l & Tech 1972,6 (1), 58 12 Westall, J C,Zachary, J L,Morel, I- M M MINEQL acomputer program for the calculahon ol chemical equd]brmm composition of aqueous systems, Water Quahty Laboratories, Ralph M Parsons Laboratory for Water Resources and Environmental Engmeerlng, Department of Civd Engmeermg, Massachusetts lnstztute ot Technology, Techntcal Note No 18, 1976 13 Sfllen, L G and Martell, A E Stabthty Constants of Metal 1on Complexes, The Chemical Sooety, London, 1964 14 Klrkner, D J , Jennmgs, A A and Thels, T L Multtsolute Mass Transport wtth Chemtcal Interactton Kmettcs, in press
A d v Water Resources, 1984, Volume 7, September
125
A numerical algonthm for predicting the migration of multiple subsurface pollutants has been developed accounting for dispersion, convection, soluble complexatlon and sohd phase accumulations (sorption) The basis of the model is a fimte element solution of the mass transport equation The essence of the algonthrn is the treatment of the sorptlon terms as implicit functions of the total soluble concentrations facilitating a general and modularized treatment of solution phase and solid phase chemistry Examples are presented dlustrating the effects of soluble complexatlon and competitive sorptlon on the transport of multIcomponent solutions
INTRODUCTION The movement of groundwater contaminants has received increased attention In recent times Both organic and inorganic pollutants from such sources as solid wastes and sludges disposed of on land, wastewaters renovated through groundwater recharge, and polluted surface waters which recharge groundwater aquifers are examples where the need for an understanding of the pnnclples governing mass transport exists Frequently complex physical, chemical, and biological interactions are Involved and must be included in the analysis of the problem. The purpose of this paper is to present a numerical algorithm to solve the equations governing the migration of subsurface pollutants accounting for dispersion, convection, soluble complexatlon and solid phase accumulations (sorptlon) Several researchers have addressed the problem of the transport of chemically reacting substances employing different methods and with varying degrees of generality For multlcomponent transport the problem is, In general, governed by a system of nonhnear partial differential equations This system of equations can be formulated explicitly If the relationships between the solid and the soluble form of the species are such that the solid phase terms can be eliminated from the transport equations Rubin and James ~ formulated the governing equations for multlcomponent one-dimensional, transport where ion exchange ~s the sorptwe mechanism The resulting equations were solved by a Galerkln finite element method (FEM) using a Crank-Nlcolson predlctor-corrector time marching scheme Valocctu et al 2 extended this technique to two-dtmenslonal domains Jennings Klrkner, and Thels 3 formulated the govermng equations for multlcomponent, one-dimensional transport Including soluble complexatIon and nonlinear adsorption These equations were solved using the FEM with a backward difference time marching scheme and a Newton-Raphson iteration within each tinle step Lln 4 formulated the problem for one-dimensional transport in a two-layered medium Including nonhnear adsorption with a Freundllch Isotherm He reduced the PDEs to ODEs by orthogonal collocation and then rotegrated usIng a fourth-order Runge-Kutta scheme Accepted April 1984 Discussion closes November 1984
120
Adv Water Resources. 1984, Volume 7, September
In all cases, above, substituting the chemical equations into the transport equations causes some or all terms to be nonhnear An alternative approach IS to omit the substitution and to treat the problem as mass transport wath a nonlinear source term This approach was taken by Trotta and Del Giudice s and Winters et al 6 in modeling coupled heat and mass transfer If the numerical scheme is such that the source term must only be evaluated but need not be known explicitly, then complicated chemical systems can be handled. In this paper the latter techmque is adopted The essence of this algorithm is the treatment of sorptive terms as tmphclt functions of the total soluble concentrations It will be shown that th~s allows a very general and modularized treatment of solution phase and sohd phase chemistry
NOTATION Bold face mlnlScules u, ~, etc represent vectors The number of components of each vector wdl be given where first used or else should be clear from the context Bold face maglscules are matnces K, M, etc The dot product of two matrices of the same order is defined A-B = Y,,IA,IBu A superposed dot indicates a spatial time derivative, u = (O/Ot)u(x, t) The gradient of a vector is the matrix Vu with components (Vu)u = Ou~/Ox/ The divergence of a vector u is defined dlv ( u ) = Z, Ou,/Ox, Lightface symbols are scalars u,, O, t, Atr etc
STATEMENT OF THE PROBLEM To help set the framework for this discussion, consider the simple three-component system with 1 2 ion exchange and soluble complexatlon, governed by the following reactions
2Cl + 62 ~ 2 ~ + c2, Ks cl +c3 ~ c l 3 , K 1
(1)
C2 + C3 ~" C23, K2
Cl and c2 are adsorbing metals and c3 is a nonadsorbmg hgand which undergoes solution phase reactions with the metals to form the soluble complexes, c13 and c23 6~ and 62 are the adsorbed forms of cl and c2 respectively For a saturated batch system a mass balance per unit volume for
0309-1708/84/030120 $2 00
© 1984 CML Pubhcanons
Multtcomponent solute transport wtth sorptton and soluble complexatton D J Ktrkner, T L Thels, A A Jennmgs each component yields after substituting the mass action equations from ( 1 )
where 0BI W 0B2 = 0B is the complete boundary and n is the outward unit normal and u(x 0)=U,o(X)
/ e l / + K I [ C 1 ] [c3] -k--~{Cl} -- [CIT]
[e~l +/,'~ [o~1 [c31 + -p0[e:IK_~( [e,] IIe,]i~ = [e~rl
(2)
[e3] +Ki[Cl] [e3] +K2 [c2] [e3] = [e3T] where [ ] = mass/volume of solution, { } = mass/mass of solids, and [c,T ] is the specified total concentration of component t Here p is the bulk density and 0 is the porosity of the domain The remmmng required equatmn is the mathematical expression of the definition of the exchange capacity, { dT}, 1 e
(¢1} + {e:) = {er}
(3)
or using the mass action equation from (1),
{if'l} + ~
\~'~!
= {CT)
(4)
Equations (2) and (4) are four equations for the four unk.nowns [cl], [c2], [ca], and [61] (they could alternatively be expressed with [c2] as a pnmary unknown) The appropriate mass achon equations then allow the soluble complex concentrations, [c13] and [c23] and the remmning sohd, {dz} to be determined For a flowing system equation (2) is simply a statement of the fact that the total concentration of each component is the sum of the species solution and sohd phase concentrations The conservation of mass equations for a flowing system must additionally be written This is facdltated by the introduction of the total soluble concentrations as variables, 1 e let (from hereon the brackets denoting concentratmn will be omitted) U 1= C1
+Kxcxe3
U2 = C2 ~/k'2c2e 3
(5)
u3= e3 + KIelC3 + K2c2e3 Now equation (2) can be rewritten I'll
+_o
(1 = ClT
.{_._PPC2 ( ( 1 ) 2 "c-~1" = c2v
0 KS
(9)
where B IS a hounded regular region Equations (5)-(9) constltute a complete statement of the problem [inphclt in the above is the local equilibrium assumption, that is the reactions, (1), are assumed fast enough compared to the transport process that they may be assumed instantaneously at equilibrium. The O~/Ot term an (7) couples the differential equations and makes them nonlinear in general This virtually rules out analytical solutions There are, however, inany approaches possible for a numencal solution Jenmngs et al 3 subsntuted equanon (5) Into (7) and also expressed (: in terms of el, cz, c3 (they considered Langmulr type adsorpnon) thus converting the governing differential equations with u as dependent variable to a differential equation with species concentration e = (c~, e2, e3) T as dependent variable This approach was also used by Valocchl et al 2 for an ion exchange system similar to that considered here only without soluble complexation The disadvantage of this technique IS that the program to solve the differential equations becomes very problem dependent This makes changing the number of components or the form of the solid phase partinonlng difficult Additionally, as the complexatlon increases (1 e as K1 and K2 become large forcing a greater percentage of the mass to exist as a nonadsorbmg complex), the solution of the differential equations becomes Increasingly difficult 3 The techmque presented here retains u as the primary dependent variable In (7) and treats ~ as an lmpllot funcnon o f u This method, as will be shown, keeps the chemical interaction calculations separate from the transport calculations and in fact allows the solution phase chemistry to be treated separately from the solid phase chemistry The method of dlscretazatlon Is the Galerkln fimte element method In space and a fimte difference method m time The resulting nonhnear algebraic equanons are solved by a Newton-Raphson procedure Other methods for discretizatlon could be employed What is emphasized here is the treatment of the chemical mteracnons due to the sorpnon term xn (7) SOLUTION ALGORITHM
0
u:
xEB
(6)
Taking the dot product of both sides of (7) with an arbitrary vector w and integrating over the domain, using the dwergence theorem and (8) yields the weak form of (7)
U 3 = C3T
The conservanon of mass equanons for saturated, lsochonc, groundwater flow with sorpnon are 7 dw (D Vu) = (Vu)v + u + (p/O) e
(7)
where u IS the vector of total soluble concentrations, 5 are the sohd concentrations, D IS the dlspersmn tensor, and v is the interstitial velomty field (dlvv = 0) For the sample problem introduced above
f w'qdA
B
i~B:
(10) where w E H ° = { w l f ( w ' w + Vw'Vw) d V < ~', w(x) = 0, x E 0BI}
B Construct finite element Interpolations for approximate solutions to (10) as follows Let
I1 = (Ul, U2, U3) T e =
f {(DVu)'Vw+w'[(Vu)v+u+(p/O)~]}dV=
(6~, d2, O)r
N'M
Statable boundary and lmtlal condltmns are reqmred with (7), for example u(x, t) = fi(x, t), x E 0BI D(Ve) n = q, x E 0B2
(8)
U ~-u h =
~
~laj(t)
1--=1 N,M
(11)
(p/O)e ~- (o/O)e~ = ~ ~lbl(t) 1=1
Adv WaterResourees, 1984, Volume 7, September
121
Multtcomponent solute transport with sorptton and soluble complexatton D J Ktrkner, T L Thets, A A Jennmgs where N is the number of nodal points and M iS the number of chemical components (in the sample problem above. M = 3) In equation (11)
¢~ : (~i. ~ ,
, ~m) T
(12)
are finite elemei~ basis functions with the following property ~]k(X~)=81t,
l=k+(t--1)M
(13)
where x~ means the tth nodal point and 6tt is the Kronecker delta Equation (13) lmphes
degrees of freedom in the standard manner Finally (17) can be wntten m vector notation as follows Ka + Ma + Mb = f
The global matrices m (19) may be assembled from element contributions m the standard manner s Equation (19) is called the semldiscrete Galerkm approximation Discretlzatton m tune further reduces the set of ordinary differential equations, (19), to algebrmc equations A vanable two-step time marching scheme is employed here s Let 1 a=--(an+l--an), At
Uk(X,) = at
(p/O)Ck(X,)=bt,
(14)
l=k +O--1)M
Thus the coefficients in the expansions (11) are nodal values of soluble concentrations and adsorbed solid concentrations This means also that only one term m each vector ~/is nonzero This is expressed mathematically as $jk(x) =
{
:~ 0, = 0.
k = l(mod M) k :/:l(modM)
(15)
N'M I [(DV~j)'V~, + ~, "(V~/v) ] dVa] • B
+
(c#,.¢])dVb]
(¢,-¢j)dVa t + ./:1
LB
* B
=
(¢,-q)dA aB:
(16) or
N.M (K,/a I +Mqa I +M,tb/) =f~,
t ~ [ t n , tn+l]
1
(20)
b=--(bn+|--bn), At
t E [ t n , tn+,]
and
a(tn÷cO = aan+l + (1 - - a ) a n
l = 1,N'M
(17)
/=i
where K,j = i [V$,'(DV~bj) + ¢,'((V¢t ) v)] dV
with 0 < a < 1 In the above
tn+ ~ = (n + a) At Evaluating (20) at tn+ ~ yields 1 l(an+' + ~tt Mbn+| = p
l(=aK
1
+--M At
1 p = f + ( I ( - - K) an + Att Mbn Various standard time marching schemes result by a suitable choice of o~ For instance a = 1 is the backward difference scheme and a -- ½ is the Crank-Nlcolson method Equation (23) is a set of nonlinear algebraic equations The simplest lteratwe scheme for solving these is to assume a guess for bn+ 1, say 0, and solve for an÷l Then solve the chemical problem at each node (equations (3), (4) and (5)) to obtain an unproved value for bn+ 1 This procedure IS then repeated until convergence Note that because of equation (15) I( and M have no off diagonal temls that represent coupling between components The only couphng that occurs is between degrees of freedom representing the same component at different nodes Thus (23)can actually be solved as M one component problems The principal &sadvantage of this approach is that It does not always converge In fact for many cases of practical importance this simple method of successwe approximations is madequate An alternative procedure is to solve (23) by a NewtonRaphson lteratwe technique 9 This can be wntten
where
],= f
(18) g = p -- I~an+l
(¢,-q) dA
~)B2 The equatJons (1 7) can be reduced to account for specified
122
Adv Water Resources, 1984, Volume 7, September
(23)
where
m+| rr/ an.l = an+| + (jm)-lgm
B
(22)
an = a(tn)
B
M,] = f (q~,'¢/) dV
(21)
b(tn+~) = abn+l + (1 - - a ) bn
It should be noted that the shape functions for the adsorbed solids do not have to be the same as for the solutes In fact since the adsorbed solids are considered implicit functions of the solutes only the inner product of the test functions and the vector of adsorbed solids has to be evaluated, the actual values of the adsorbed solids are never required There are two reasons for adopting the interpolation scheme of equation (11), however First the amounts of the adsorbed solids are of intrinsic interest Second a substantial amount of computer time can be spent evaluating the inner product mentioned above, since it will be shown to involve chemical speciatlon When the solids are interpolated as m (1 I) the number of chemical calculations is kept moderate whde also reducing the inner product calculations since the M matrix (see below) only needs to be computed once Now using (1 1) m (10) and taking the test functions w to be ~ (t = 1, N . M ) yields
/=1
(19)
1 ~tt Mbn+|
1 J : r, + - -
At
MJ
J = (ab/Oa) [tn+,~
(24)
Multwomponent solute transport with sorptlon and soluble eomplexatton D J Ktrkner T L Thets, A A Jennmgs and the superscript m ln&cates the iteration number The principal difficulty w~th this method hes m the calculation of the Jacoblan matrix, J, since b is given as a function of a only implicitly Note that by virtue of equation (15), J has a block dmgonal form (see Fig 1) where each submatnx is itself a Jacobmn matrix of the form Jt = (Oe/OU) at (xt, tn.~),
l = 1,N
(25)
and each J~ is of dimension M x M Since e is an implicit funcnon of u, one approach for evaluating Jt is to calculate the derivatives numerically lo Oat
at(Ut,
, U1 + 6,
,Um)
0u /
-
-
Ct(U 1,
, U 1,
, UM )
Thus e = goh-l(u)
(30)
(ae/0e) Oe/Ou)
(3 l)
and then from the chain rule
~e/0u
=
Since equanon (28) is often solved numerically m chemical equ~hbrlum programs by Newton's method, the Jacoblan matrix 0h/0c is available lz From equanon (28)
ac/0u =
(0h/be)-'
(32)
For many exchange expressions (such as the mulncomponent Langmulr isotherm), the matrix 0e/0e may be evaluated explicitly Thus equation (31) may be written
e
Oe/au = (~e/Oe) (0h/0e) -t
where e is some small parameter To evaluate each term on the right requires (a) Solution phase speclanon For the sample chemistry introduced previously th~s means solving equation (5) for ct. c:, c3 gwen u (b) Sohd accumulation Given c~, c2, cs, use equation (4) and the mass action equation for (1) to calculate 61 and (2 The key aspect of the above algorithm ~s that the chemical computations may be modularized, 1 e kept ennrely separate from transport calculanons In fact solution phase computations are distract from sohd phase computations thus facilitating the study of different sorptlve mechanisms Some degree of generahty may be obtained ff the solunon phase specmnon calculations (part (a) above) are performed by one of the available general-purpose chemical equthbna programs (e g Morel and Morgan H) It must be noted that when a fairly comphcated chemical module is employed the computanonal rime reqmred to form the Jacobmn matrices, Jt may become rather burdensome Use of a modified Newton iteration, I e use jo instead of j m in equation (24). reduces computatmnal nine considerably and has been found to work well An alternative method, however, may be employed to calculate Jt by using the reformation avadable m the solunon phase chenucal module First, to generahze equahon (5) define the vector of species concentranons as
c = (~ l,
= g(h-l{u))
, CM)T
(27)
Then equanon (5) can be written u = h(e)
(28)
where h Is m general a nonhnear operator Equation (4) and the mass action equatnm tar (1)~ can be generahzed to = g(c)
Equatmn (33) may be used in lieu of equation (26) Equation (26) requires M + 1 calls to the chemical module of the program to evaluate J . whereas equatmn (33) only requires one call However, (33) requires the inverse of 0h/0c The reverse of 0h/0c nonetheless is easier to compute than speclatmg M times Also note that equations (33) or (26) are best computed globally and thus the global matrices, K and M need only be computed once Results of a sample problem wtll now be presented to demonstrate the type of problem that can be addressed with th~s algorithm
EXAMPLE
In ref 1 Rubm and James studied several systems where ion-exchange was the sorptlve mechanism One particular example involved the leaching of a sodlc soft governed by a 1 2 binary exchange reaction The transport is one-dimensional with a mixed or tlurd type boundary condmon at the source and a zero gradient at the downstream boundary (to s~mulate a semHnfimte region) This hypothencal problem wdl be reexamined here and the effect of the addmon of a nonadsorbmg hgand (SO~) which forms soluble complexes with the adsorbing metals will be stu&ed It should be emphasized that this example is not offered as a thorough analysis of the Ca/Na/SO4 transport system, but rather to dlustrate the potential of the method presented hereto to study such systems Thus the governing reactions are equation (1) with cl = Na +, ~2 = Ca m and c3 = SO4 The mass balance and boundary conditions for thts problem are 02tl~
OU t
D--=v Oxz
Ox
+u,+(p/O)6,,
t=1,3(¢3=0)
(34)
I = 1,3
(35)
(29)
( --D'~--x Ou, + vu,} = v~,, "lx=O Out
DEGREE OF FREEDOM NODE
1
I I
2
f I
. . . . . . .
(33)
i I
x=r=0'
1=1.3
(36)
N
The m~ttal condlhons are
co--iA = Jl,,,oEJ21-i,.., o., , liON] Ftgure 1 Contrtbutton to the Jacobtan matrix oJ the sorptton terms
ut(x,O)=uto,
1=1,3
(37)
Numerical values for all the physical ~.onstants, boundary condmons and chenncal data are given m Table 1 Rubm and Janles considered three cases for this problem The firs! one (ref 1 case H) had a high calcmm lnfluent and no sulfate (U3o = u3 = 0 0) The second case (ret 1 case L) was identical in all respects to the first case except the calcmm mfluenl concentration was low The
Adv Water Resources, 1984, Volume 7, September
123
M u l t w o m p o n e n t solute transport with sorptton and soluble complexatton D J Ktrkner, T L Thezs, A A Jenmngs Table 1
Datafor sample problems
160
ut
ut0
Component
1 2 3
0 0020 0 0072 0 100
0 0280 0 0020 00
Na Ca SO,
1
I
I
I
140
Boundary and amtml conditions (meq/cm 3) l
I
120 100
u,
80
/ /A / /...'o.,
G, 6 0 40
Physical data Symbol
Value
Umts
0
0 54 1.20 25 36 2 536
g/cm 3 cm2/day cm/day
....
P'~ '--~-
D v Chemical data Symbol
Value
K, K: KS
2 22 79 2 0 078
o
1' 0
I
I
40
80
1
I
I
120 160 D E P T H , cm
7'00
240
Figure 2 Normahzed sodium concentration versus depth Ul - total soluble sodium concentration, u2 - total soluble sulfate concentration
Umts cm3/meq * cm3/meq * cm 3/gt
I0 U z w, t h ~ 3
~
08 '~
* Equfl]brmm constants trom ref 13, corrected to an average constant xomc strength of 0 05 using the Davies equation From ref 1
\
third case was similar to the second case except gypsum was present in the soft It was assumed that the gypsum was present in a sufficient amount so that the calcium-sulfate solubility product was always sausfied This allowed the number o f operational equations to be reduced to two Herein the first case is reexamined, but w~th sulfate added to the mfluent source Sulfate has been added at a concentration that will allow for complexaUon, but not precxpltatlon o f calcium sulfate Thus there are still three operational equaUons (equation (34)) This example serves to demonstrate the application of the algorithm presented m this paper to a problem with soluble complexation and adsorption Linear basis funcuons are employed for the fimte element Interpolation o f solutes and solids, and a backward difference time marching scheme is used Figure 2 shows normahzed concentration profiles for total soluble sodium at 14 6 and 43 8 days The sohd lines are results with no sulfate mfluent, agreeing with Rubln and James' case H The dashed lines are the total soluble sodmm profiles when the sulfate mfluent is 0 1 0 m e q / c m 3 The presence o f the sulfate has a small effect on the sodium profiles Figure 3 shows normahzed concentration profiles for total soluble calcium w~th and w~thout sulfate present and for total soluble sulfate The presence o f the sulfate causes a calcmm pulse to propagate thiough the domain along with the sulfate front This effect can be explained as follows Figure 4 shows several batch isotherms for this chemical system, 1 e adsorbed c~ (or c:) versus c2 for fixed c~ F o r the simulation presented above, c~/clo is always less than 14 and thus the soil has a strong affinity for c2 except at very low concentrations o f c2 (Ca 2+) Given the original pore water concentrations the behavior o f the profiles without sulfate is to be expected Ion exchange only affects transport in the mmal portion o f the domain
Adv Water Resources, 1984, Volume 7, September
=00
u z w l l h "~5 = 0 ' 0 .
~
~
'
-
U3
-
\
U, 04
124
0 ~
20 o
p
~u 3 " :
.~.~"
ok]:~--;~~--0
40
,\''"r'--~
,
80
120
160
DEPTH,
200
240
cm
Figure 3 Normahzed sulfate and calcium concentration versus depth u ~ - total soluble sodium concentration, u2 - total soluble calcium concentration, u3 - total soluble sullate concentration
I0-
)
CI /C)o= 5
i
i
I
0
.
08
02
06
04
Ez
T,
CT 04
06
02
08
0
0
I
02
I
04
i
06 C2/C20
I
08
TT
I0 0
Figure 4 Sorbed calcium versus soluble calcium concentration cl - [Na÷], (1 - sorbed sodium, c2 -[Ca++], c2 sorbed calcium, ET - cation exchange capacity
M u l t t c o m p o n e n t solute transport wtth sorptton and soluble complexanon D J Kwkner, T L Thezs, A A Jennmgs
0 8
principal advantage o f this approach is the modularlzanon o f the chemical calculanons Solunon phase calculanons are also performed separately from solid phase calculations Add]nonal one- and two-dimensional mult]component apphcat~ons with complex chemical interactions are m preparation It should also be mentioned that tMs approach readily extends to the problem o f mass transport with sorpnon governed by chemical kmet]cs 14
8 Ooys 146 Doy s - ' ~ cz I u z
c~., un
~ C l /
UI
04
ACKNOWLEDGEMENT 02
0
I
40
80
I
I
120 160 DEPTH, cm
I
I
200
240
2BO
Ftgure 5 Fracture o f each c o m p o n e n t as free ton versus depth c] - [Na÷], u i - total soluble sodtum concentratton, c2 - [Ca*+], u 2 - total soluble calcium concentratton
V,qth sulfate present the calcmm (uz) in solunon does not all exist as Ca 2÷ (c2) and the sodium ( u 0 m solution does not exist only as Na ÷ (cO because o f c o m p l e x a n o n With sulfate present the calcmm (sodium) ex]sts as the exchanging free ion, Ca 2÷ (Na*) and as the complex CaSO~(aq) (NaSO~) Figure 5 shows a plot o f the free ion concentration o f each c o m p o n e n t as a fraction o f the total soluble concentration versus depth at 14 6 and 43 8 days Notice that not much sodmm is complexed However, m regions o f Mgh sulfate concentranon (compare with F~g 3) a large percentage of the calcmm is complexed Recalling the Isotherms in Fig 4 it is expected that adsorbed sodmm decreases while adsorbed calcmm increases Since this ~s 1 2 exchange, one should expect a greater increase m the total soluble sodmm than the corresponding decrease m total soluble calcmm This is indeed the effect displayed m Figs 2 and 3 In the original smmlanon w~thout sulfate, ion exchange only affects the first 10 cm o f the domam With sulfate added, a slgmficant degree o f soluble complexatlon occurs (especially CaSO4°) which affects ion exchange as far into the domain as the active sulfate front has penetrated CONCLUSIONS In this paper, an algorithm has been presented for solving the mass transport equanons with sorpnon and soluble c o m p l e x a n o n The procedure employs a Galerkm finite element spanal d~scretlzanon and a finite difference time marching scheme S o r p n o n terms are treated as lmphclt functions o f the total soluble concentranon variables The
The work presented m this paper was supported by grant DE-AC02-79EV 10253, from the Ecological Research Dwlslon, Office o f Energy Research, US Department of Energy to the University o f Notre Dame
REFERENCES 1 Rubln, J and James, R V Dlsperslon-attected transport of reacting solutes in saturated porous medm Galerkm method apphed to equxhbrmm controlled exchange m umdlrechonal steady water flow, WaterResour Res 1973, 9 (5), 1332 2 Valocchl A J , Street, R L and Roberts, P V Transport ot ion-exchanging solutes m groundwater chromatographic theory and held simulation, Water Resour Res 1981, 17 (5), 1517 3 Jennmgs, A A , Ktrkner, D J and Thels, T L Multi-component equfllbrmm chemistry m groundwater quality models, Water Resour Res 1982, 18 (4), 1089 4 Lln, S H Nonlinear adsorption m layered porousmed]a flow, J Hydrauhcs Dlv, ASCE 1977, 103 (HY9), Proc Paper 13192, 951 5 Trotta, A and Del Gmdlce, S A finite element solution of coupled heat and mass transfer with chemical reaction, Chem Eng Sct 1980, 35, 709 6 Winters, K H , Rae, J , Jackson, C P and Chffe, K A The finite element method for laminar flow with chemical react]on, Int J Num. Meth k.ngng 1981, 17, 239 7 Bear, J Dynamws o f Fluuts m Porous Medta, American Elsevier Co , New York, NY, 1972 8 Zlenk]ew]cz, O C The Fmtte Element Method, McGraw-Hall, 1977 9 Ortega, J M and Rhemboldt, W C ]terattve Solutton of Non hnear Equattons zn Several Vartables, Academic Press, New York, 1970 10 Oden, J T Fmtte Elements of Nonhnear Contmua McGrawHdl, 1972 11 Morel, F and Morgan, J J A numerical method for computing equlhbna m aqueous chemical systems, finvzron S c l & Tech 1972,6 (1), 58 12 Westall, J C,Zachary, J L,Morel, I- M M MINEQL acomputer program for the calculahon ol chemical equd]brmm composition of aqueous systems, Water Quahty Laboratories, Ralph M Parsons Laboratory for Water Resources and Environmental Engmeerlng, Department of Civd Engmeermg, Massachusetts lnstztute ot Technology, Techntcal Note No 18, 1976 13 Sfllen, L G and Martell, A E Stabthty Constants of Metal 1on Complexes, The Chemical Sooety, London, 1964 14 Klrkner, D J , Jennmgs, A A and Thels, T L Multtsolute Mass Transport wtth Chemtcal Interactton Kmettcs, in press
A d v Water Resources, 1984, Volume 7, September
125