A numerical model of coupled fluid flow and reactive geochemical transport

607

A numerical model of coupled fluid flow and reactive geochemical transport Gour-Tsyh(George) Yeh a and Yilin Fang b aDepartment of Civil and Environmental Engineering University of Central Florida, Orlando, FL 32816-2450, USA bDepartment of Civil and Environmental Engineering The Pennsylvania State University, University Park, PA 16802-1408, USA Experiments have shown that geochemical reactions have profound effects on the mechanics of fluid flows. Changes in hydraulic conductivity and fluid pore velocity were observed in the unsaturated flow experiments. These changes were attributed to secondary phase precipitation, which occurred under reduced saturation conditions. Clearly the depletion of water by precipitation would manifest itself if the water is not supplied via hydrologic transport. It is not difficult to see that the heterogeneous reactions including precipitation can plug pores reducing both the porosity and hydraulic conductivity, which in turn would hinder the transport of water into the affected regions. As a result, moisture content would be reduced and further precipitation would occur. The onset of further precipitation would remove water from the liquid phase. These complicated interactions between flow dynamics and geochemical reactions must be included in modelling reactive transport, especially when the degree of saturation is low. This paper communicates the development of a mechanistic-based numerical model for simulation of coupled fluid flow and reactive geochemical transport, including both fast and slow reactions, in variably saturated media. A simple example is used to illustrate the interplay and effects of advective/dispersive transport, geochemical reactions, and fluid flows. With the ability to simulate coupled fluid flow and reactive transport, an effective strategy for investigating contaminant transport and for conducting performance assessment of various waste management and remediation technologies can be developed. Keywords: reactive transport, numerical models, groundwater flow, contaminant transport. 1. I N T R O D U C T I O N The couplings among chemical reactions, advective and diffusive transport in fractured media or soils, and changes in hydraulic properties due to heterogeneous reactions within fractures and in rock matrix are important for both waste disposal and remediation of contaminated sites. Experiments have shown that geochemical reactions have profound effects on the mechanics of fluid flows. Rapid, dramatic changes in hydraulic conductivity and fluid pore velocity were observed in the unsaturated flow-through testing in columns

608 filled with glass beads [1]. These changes were attributed to secondary phase precipitation, which occurred under reduced saturation conditions. Clearly the removal of water by precipitation would manifest itself if the water is not supplied via hydrologic transport. It is not difficult to see that the heterogeneous reactions including precipitation can plug pores reducing both the porosity and hydraulic conductivity, which in turn would hinder the transport of water into the affected regions. As a result, moisture content would be reduced and further precipitation would occur. The onset of further precipitation would deplete water from the liquid phase. These complicated interactions between flow dynamics and geochemical reactions must be included in modeling reactive transport, especially when the moisture content is low. The development of mechanistically-based reactive chemical transport models has exploded in the last decade [2,3]. Some models couple transport with equilibrium chemistry, others couple transport with kinetic geochemistry for specific geochemical processes like precipitation-dissolution or adsorption. This paper describes the development and application of a mechanistically-based numerical model for simulation of coupled fluid flow and reactive chemical transport in variably saturated media. The model include a full suite of kinetic and equilibrium geochemical reactions to address the complexity of geochemical processes in the subsurface environment. 2. G O V E R N I N G E Q U A T I O N S The governing equations of the model include three modules: a flow equation, a set of hydrologic transport equations, and a set of geochemical equations. A modified Richard equation describes density-dependent fluid flow in variably-saturated media. It can be derived based on continuity of fluid, continuity of solid, motion of fluid (Darcy's law), consolidation of the media, and compressibility of water as well as depletion of water via precipitation reactions. The governing equations for a reactive system are mass balance equations fore every chemical species, which can be reduced to a set of transport equations for kinetic variables and mass action equations [4]. For simplification, we will present the governing equations for fluid flow and reactive chemical transport without including source/sinks; detailed descriptions of these equations can be found elsewhere [4,5].

2.1. Governing Equation for Fluid Flow The derivation of the governing equations for flow through saturated-unsaturated media can be derived based on (1) continuity of fluid, (2) continuity of solid, (3) motion of fluid (Darcy's law), (4) consolidation of the media, (5) compressibility of water, and (6) removal of water due to geochemical reactions as [5] .Po Ot.

- V. . V -.

P ;. V

. K

.

V h. + Vz. ; F

a' ne Oe + ~'Oe + ne~-~ dS

(1)

where 0e is the effective moisture content (L3/L3), h is the pressure head (L), t is the time (T), z is the potential head (L), W is the removal rate of water due to heterogeneous chemical reactions such as precipitation (M3/L3/T), po is the referenced fluid density at zero chemical concentration (M/L3),p is the fluid density (M/L3), a' is the modified compressibility of the soil matrix (l/L),/3' is the modified compressibility of the liquid (l/L), ne is the effective porosity (L3/L3), and S = degree of saturation.

609

2.2. Governing Equations for Reactive Chemical Transport Starting from the mass balance equations for the system of M chemical species and N geochemical reactions, we can decompose the reaction matrix to yield [4] 00E~ 00E? 00El' = L(E~)+ DkkRk, i e MKI, k e NKI 0t + 0t + 0---p-

(2)

K~ =

(3)

(Ai) "'k /

(Ai) tqk , k e NE

where 0 is the moisture content; E~, E~, and E~ are the total dissolved, sorbed, and precipitated concentrations of the i-th mobile kinetic variables, respectively; L is the hydrologic transport operator; Dkk is the k-th diagonal entry of the decomposed reaction matrix; Rk is the reaction rate of the k-th kinetic reaction; MKI is the number of kinetic variables; NKI is the number of linearly independent kinetic reactions; K~, is the equilibrium constant of the k-th fast reaction; Ai is the activity of the i-th chemical species; Vik is the reaction stoichiometry of the i-th species in the k-th reaction associated with the products; Pik is the reaction stoichiometry of the i-th species in the k-th reaction associated with the reactants; and NE is the number of equilibrium reactions (fast reactions). It should be noted that when Dkk is equal to zero, the i-th kinetic variable Ei becomes a chemical component. It is further noted that MK! -- NKI -+- Nc =- M - NE, Nc being the number of components. 2.3. Coupling between Fluid Flows and Reactive Chemical Transport The density-dependence of fluid flows are described via the following constitutive relationships MKI

MKI

P = Po + E E~Mi, # = #o + E aiE~'Mi i=l

(4)

i=l

where p and # are fluid density and fluid dynamic viscosity; Po and #o are fluid density and fluid dynamic viscosity in pure water; Mi is the molecular weight of the i-th aqueous kinetic variable; and ai is a weighting parameter for the i-th aqueous kinetic variable. The interactions between reactive transport and fluid flows are described by the following dependencies of moisture content (0), saturated hydraulic conductivity tensor (Kso), and dispersion coefficient tensor (D) on concentrations of minerals and the dependency of water depletion rate on the rate of production of all precipitated species

SOlo

0 = 1 + S~p

and K . o = K . o o (

1 )n 1 + S~p

D = Do [0m-l(1- ~p)m] and Mw = ~

aiwr~

(5)

(6)

iEMp

where 0.~ois the effective saturated moisture content without the effect of solid species; S is the degree of saturation of water which is a function of the pressure head; ~p - ~ piVi, where Pi is the precipitated concentration of the i-th mineral [mole/dm 3 of water] and Vi

610 is the molar volume of the i-th mineral [dm 3 of solid/mole]; Ksoo is the reference saturated hydraulic conductivity tensor (ignoring the effect of density and solid concentration); n is the fractal for estimating hydraulic conductivity based on particle size and packing structure; Do is dispersion coefficient tensor computed without the effect of solid species; m is the "cementation exponent" [6]; Mp is the number of precipitated species; aiw is the stoichiometric coefficient of water in the i-th precipitated species; and r~ is the production rate of the i-th precipitated species. 3. N U M E R I C A L

SOLUTIONS

The flow equation is discretized with the Galerkin finite element method. The transport equations are discretized with either conventional finite element methods or a hybrid Lagrangian/Eulerian finite element method. Nonlinearities in the fluid flow, reactive transport, and their interaction are solved with using Picard linearization, and any nonlinearity in geochemistry are solved with Newton-Raphson linearization. The two subsystems of equations, one describing the hydrologic transport of mobile kinetic variables and the other describing the mass balance of immobile kinetic variables and fast equilibrium reactions, are solved iteratively. Three different approaches are to reach a convergent solution. The first approach is complete iteration between two subsystems, the second approach is the use of operator splitting, and the third approach is the employment of the predictor-corrector method. For each time step, a sequential iteration approach is also applied to solve coupled fluid flow and reactive chemical transport equations until the pressure head and total analytical concentrations of mobile variables meet prescribed convergent criteria. 4. E X A M P L E This simple problem simulates the effects of matrix diffusion and precipitation reaction induced liquid removal on solute and reactive transport through fractured media. The model was based on a longitudinal cross-section of a 0.76 meter high column with a 120 #m wide fracture in the center, running the length of the column. Due to symmetry, one-half of the cross-section is used; the resulting domain (Fig. 1) is 7.6 dm high by 0.325 dm wide with a 6 x 10-4 dm wide fracture on the left hand side. The fractured media is simulated by using specific elements having a moisture content, 8, of 1.0. The matrix elements have a moisture content (without the presence of precipitated species), tortuosity, and bulk density of 0.1, 0.1, and 1.7 kg/dm3, respectively. The referenced saturated hydraulic conductivity Ksoo in the z-direction for the matrix is Kz -- 2.16 x 10 -6 dm/hr; in the fracture, Kz = 237.0 dm/hr for material 2, K~ = 222.0 d m / h r for material 3, and K~ = 90.0 dm/hr for material 4. The diffusion coefficient, D, is 3.6 • 10 -4 dm2/hr in the fracture cells and 3.6 • 10-~ dm2/hr in the matrix cells. Longitudinal dispersivity, OlL for all elements is set to 7.6 dm while lateral dispersivity, aT, is set to 0.0 dm for all elements. For flow simulations, both top and bottom surfaces are held at a Dirichlet head of 7.615 dm and 7.6 dm, respectively. The matrix is homogeneous and is designated as Material No. 1. For reactive transport simulations, the upper boundary of the domain is assumed to maintain a constant concentration, Dirichlet boundary. Neumann boundaries are used for the other three boundaries with 0 C / 0 x = 0 for the two vertical boundaries

611

Fracture

Flow: Dirkhlct Head = 7.615 dm T r a n s p ~ Didchlet C = C o

~al

:~

22 21

ua~

352 336374 357

63 u

2O 19

'

356

3,~

~ 6.0E,4din ~ 6J~4dm ~ 13.E-3dm

15 mm of r mcnmmt|m size dov~ld

3

l~t~ill 1

7.6 dm

4

378 3

339

2

338

3 0.325 dm

z !i

2 1

Flow: Dirichlet Head = 7.6 dm Tnuupm: Ymablr B.C.

Figure 1. Problem Definition

0.8 dm

377 0.8 dm

376 I

X

396 395 394 393

22

23 ~".~-~ ,~ 1.2E-4din

43

45 ~ =~ ~ 2.4F~4din

M

316 =

IIr iac~,min I in size to the rilht

331

337 353

0.8din

375

9 -" 0.075 dm 0.1 dm

Figure 2. Finite Element Discretization

and 0C/0z - 0 for the bottom boundary. The grid is made up of 21 rows by 17 columns of rectangular elements, with the three left-hand columns of elements representing the fracture and the remaining elements representing the matrix (Fig. 2). The fracture element widths (horizontal direction) are 1.2 x 10 -4 dm in the first column and 2.4 x 10 -4 dm in the second and third columns. The matrix element widths gradually increased from being equal to the fracture element width, 2.4 x 10 -4 dm, at the fracture side to 0.1 dm at the opposite side. The heights (vertical direction) of rows of elements gradually increase from 6 x 10 -4 dm to 0.2 dm, from the top to the bottom. Three cases are simulated in this example: in the first case (Case A), the flow and reactive chemical transport are decoupled; in the second case (Case B), the flow and reactive transport are coupled through the effect of precipitation on hydraulic conductivity, porosity, diffusion, and dispersion; and in the third case (Case C), an additional interaction between flow and transport is to assume that precipitation reactions will deplete stoichiometriclly water. To reflect the effect of mineral precipitation on flow, the fractal exponent for adjusting the hydraulic conductivity (Eq. 5) is set equal to 2.5. The cementation exponent in Archie's law is set equal to 2.0 to adjust the dispersion coefficient (Eq. 6). Molar volumes (Vi) of 0.03 dm3/mol and 0.037 dm3/mol are assumed for portlandite, Ca(OH)2 (s), and calcite, CaCO3 (s), respectively. To reflect the effect of water depletion via precipitation reaction, a stoichiometric coefficient of water in portlandite and calcite are hypothetically set to 6. In Case A, although the effects of precipitation/dissolution on flow, conductivity and dispersion are not considered, solute transport and chemical reaction still need to be coupled because heterogeneous reactions (precipitation/dissolution reactions) are involved in the system. A time step size of 0.24 hour was used; a total of six days (144 hours) is simulated using 600 time steps. The chemistry of the system is relatively simple: the chemical components of the systems are Ca 2+, CO 2-, H +, and a conservative tracer. Calcite is present in the matrix

612 (8.642 • 10-5 M) in equilibrium with the pore water and fracture water which contain 10 -3 M Ca 2+ and CO ]- (total concentrations). During the simulation, a solution 10xricher in both total calcium and total carbonate is injected in the media via the top boundary. The complexed species include OH-, CaCO3, CaliCO +, CaOH +, HC03, and H2CO3. Calcite (CaCO3(s)) is included as solid species subject to kinetic precipitation and dissolution. The reaction stoichiometries, equilibrium constants and rate constants can be found elsewhere [7]. The results of the simulation are summarized by the profiles of calcite shown in Figs, 3, 4, and 5, respectively, for Cases A, B, and C, respectively. In Cases B and C, the diffusion transport into the matrix is diminished due to the calcite precipitation which impedes the diffusion of calcium and carbonate, resulting in less calcite precipitations. This is confirmed by comparing Figs. 4 (for Case B) and 5 (for Case C) with Fig. 3 (for Case A). Indeed the concentration profiles of calcite at z = 0.0024 dm for Cases B and C in the matrix zone are much smaller than those in Case A. On the other hand, since less amount of contaminant can diffuse into the matrix zone, a larger amount of contaminant will tarry in the fracture zone, in which the advective transport dominates. This can be found as the concentrations of the conservative tracer in the fracture zone for Cases B ( middle frame in Fig. 6) and C (right frame in Fig. 6) are higher than those in Case A (left frame in Fig. 6). Examination of the simulated velocity field (not shown) shows that the change of flow velocity is very small because the calcite precipitation in the fracture zone fills only a small portion of the fracture. In the matrix, the velocity is so small that diffusion transport of chemicals is dominant for all three cases. In other words, there is little hydrologic transport of water into the matrix. This should explain the difference in concentration profiles of calcite for Cases B and C: the concentration of calcite in the matrix for the later is larger than those for the former (compare Figs 4 and 5). This is so because precipitation reactions deplete water, which decreases the moisture content. As the water is not replenished, the moisture decreases and the concentrations of Ca 2+ and COl- increase to such a degree that they are saturated with respect to the calcite. This naturally induces more calcite precipitation. It can be expected that these processes can manifest themselves and may eventually plug pores, which would make the matrix impermeable and immobilize the contaminant.

5. C O N C L U S I O N This paper communicates the development of a mechanistically-based coupled fluid flow and reactive geochemical transport model including both fast and slow reactions in variably saturated media. The theoretical bases and numerical techniques were summarized and a simple problem was presented. The problem examined the effects of precipitation on fluid flow and matrix diffusion in fractured porous media. The precipitation reactions in the matrix hinder the diffusion of contaminant into the matrix. This decreases matrix retardation, which causes fast transport in fractures. The depletion of water in the matrix via precipitation reactions promote further precipitation when water is not replenished via hydrologic transport. These processes many manifest themselves and immobilize contaminants.

613

A

omTs

Dk*-nee h~m top of eotumm ~. o~o=4ck,,

o.oo~m

o

o.omtem

I o.ee5

o.otesm

i

J ....

o.~,. . . .

o'.,' ' ' ' o . ~ ' ' ' ' o~, . . . . o ~ . . . . ObBnel lrom Center M Fracture (din)

o'~

....

....

o.~ . . . .

o'.,' ' ' ' o . ' , ' ~ ' 'o~ . . . . o~ . . . . DIMancefrom Cen~r of Fracture (den)

o'~ . . . .

Figure 3. Without considering precipi- Figure 4. Considering precipitation eftation effect fect

o.ol

Ol,,kmtm from top ol ookim~ ~ o.oo~dm o.olram

o~o7= o.oo6

o=~-_~ ~ .~. .~. .~. .o~oozs o

. . . .

,

o.o6

. . . .

.~. ,

o.1

o . . . .

,

o.1 s

n . . . .

i

o~

. . . .

[] ,

o~

. . . .

,

0.3

. . . .

,

o.~s

olMmncmfrom centre' of Fracture (clot)

Figure 5. Considering water removing effect due to precipitation.

7.6

Z

7.6

Z

7.6[ Z

5.7

5.7

5.7

3.8

3.8

3.8

1.9

1.! Tracet (") 0.95

1.9 Tracer

8

Tracer (M)

(M)

O.~J5

Figure 6. Tracer concentration profile at x=0.0dm

8.9

o.~s

o~=

614

ACKNOWLEDGEMENTS Research is supported, in part, by National Science Foundation under Award No. EAR0196048 with University of Central Florida., and, in part, by U.S. Environmental Protection Agency under grant No. R-82795602 with University of Central Florida. REFERENCES

1. B.P. McGrail, C.W. Lindenmeier, P.F. Martin, G.R. Holdren, and G.W Gee, 1996. "The Pressurized Unsaturated Flow (PUF) Test: A New Method for EngineeredBarrier Materials Evaluation." In Proc. American Ceramic Society Conference held in Indianapolis, Indiana, April 14-18, 1996. 2. G.T. Yeh and V.S. Tripathi, 1989. A Critical Evaluation of Recent Developments in Hydrogeochemical Transport Models of Reactive Multichemical Components. Water Resources Research, 25 (1): 93-108. 3. C.I. Steefel and P. Van Cappellen (Guest Editors), 1998. Special Issue: Reactive Transport Modeling of Natural Systems. Journal of Hydrology, 209 (1-4): 1-388. 4. G.T. Yeh, 2000. Computational Subsurface Hydrology Reactions, Transport, and Fate of Chemicals and Microbes. Kluwer Academic Publishers. Boston, Massachusetts. 318 pp. 5. H.P. Cheng and G.T. Yeh, 1998. Development of a three-dimensional model of subsurface flow, heat transfer, and reactive chemical transport: 3DHYDROGEOCHEM. J. Contaminant Hydrology, 34:47-83. 6. F.A.L. Dullien, 1979. Porous Media. Academic Press. 574 pp. 7. G.T. Yeh and Y.L. Fang, 2002. A Coupled Fluid Flow and Reactive HYDROGEOCHEMical Transport through Saturated-Unsaturated Media- Version 3.0. Technical Report. Department of Civil and Environmental Engineering, University of Central Florida, Orlando, FL 32816.