United States
Environmental Protection
Agency
National Risk Management
Research Laboratory
Ada, OK 74820
Research and Development
EPA/600/SR-98/159
July 1999
&EPA Project Summary
3DHYDROGEOCHEM:
A 3-Dimensional Model of
Density-Dependent Subsurface
Flow and Thermal
Multispecies-Multicomponent
HYDROGEOCHEMical Transport
Gour-Tsyh (George) Yeh and Hwai-Ping (Pearce) Cheng
A three-dimensional finite-element
numerical model designed to simulate
reactive chemical transport in
subsurface systems with temperature
effect taken into account is presented.
The three-dimensional model is
developed to provide (1) a tool of
application, with which one is able to
deal with a variety of real-world problems,
(2) a tool of education, with which one
can study how a factor would affect the
whole system, and (3) a substructure,
which one could modify to handle
specific problems. The hydrological
environment to which the model can be
applied is a heterogeneous, anisotropic,
saturated-unsaturated porous medium
under either transient-state or steady-
state flow conditions. In addition, the
temperature within the system of interest
can be both time- and location-
dependent. For steady-state
simulations, strong coupling among
subsurface flow, reactive chemical
transport, and heat transfer is used in
the model. For transient-state
simulations, weak coupling is used, but
a density effect is still considered in
computations. Both the strong and the
weak couplings are pictured and
explained in the original report. The
model employs chemical equilibrium to
describe the relationship among
chemicals. The chemical reactions
included in the model are aqueous
complexation, multi-site adsorption/
desorption, multi-site ion-exchange,
precipitation/dissolution, redox, and
acid-base reactions. To discretize the
domain of interest appropriately, the
element used in the model can be a
hexahedral, a triangular prism, or a
tetrahedral element. To extend its
applicability to more real-world
problems, two approaches are presented
for the reactive chemical transport
module. The first approach uses the
pore velocity and dispersion coefficient
to handle advection and dispersion,
respectively, for aqueous components,
whereas the second approach employs
the retarded pore velocity and the
retarded dispersion coefficient. Both
approaches are designed to handle the
problems with dominating precipitated
species involved. In the original report,
the governing equations of subsurface
flow, reactive chemical transport,
chemical equilibrium, and heat transfer
are stated and/or derived, followed by
the numerical strategies for solving the
governing equations. In addition, twenty-
five designed examples used to verify
the model were described. Four
application examples, including a three-
dimensional subsurfacef low example, a
three-dimensional reactive chemical
transport example, a three-dimensional
heat transfer example, and a three-
dimensional coupled flow-transport-
transfer example, are presented to
demonstrate the capability of the model.
-------
The modeling experiencefromthis study
is summarized.
This Project Summary was developed
by EPA's National Risk Management
Research Laboratory's Subsurface
Protection and Remediation Division,
Ada, OK, to announce key findings of
the research project that is fully
documented in a separate report of the
same title (see Project Report ordering
information at back).
INTRODUCTION
3DHYDROGEOCHEM (a3-Dimensional
coupled model of HYDROIogic transport,
GEOCHEMical equilibria, subsurface flow,
and heat transfer in reactive multi-
component/multi-species systems) has
been developed mainly to study reactive
chemical transport in subsurface systems.
It is composed of two basic modules: flow
and transport. In the flow module, the
Galerkin finite element method is used to
discretize the modified Richards' equation
for simulating density-dependent
subsurface flow. In the transport module,
both reactive chemical transport and heat
transfer are simulated. Their governing
equations are discretized by using either
the conventional Eulerian finite element
method or the hybrid Lagrangian-Eulerian
finite element method. In reactive chemical
transport, solute transport is coupled with
chemical equilibrium to consider the
interaction among chemical species as they
are transported by subsurface flow.
3DHYDROGEOCHEM has the following
features which make it flexible and versatile
in modeling a wide range of real-world
problems: (1) Irregularly-shaped domains
can be dealt with; (2) Both heterogeneous
and anisotropic media, as many as desired,
can betaken into account; (3) Both steady-
state and transient-state simulations can
be processed; (4) The initial conditions can
be either prescribed or obtained by
simulating a steady-state system under
consideration; (5) Both distributed and point
sources/sinks can be considered spatially-
and/or temporally-dependent; (6) All types
of boundary conditions can be considered
spatially- and/or temporally-dependent; (7)
The appropriate variable boundary
conditions are determined automatically;
(8) The "in-element" particle tracking
technique is used to accurately and
efficiently perform particle tracking; (9)
Either the conventional Eulerian or the
Lagrangian-Eulerian finite element method
can be used to solve transport equations;
(10) A numerical scheme of capturing peaks/
valleys is employed to increase the
computational accuracy of reactive
chemicaltransport; (11) Multi-adsorbing site
and/or multi-ion exchange site options are
available to accountfor chemical reactions;
(12) Many options are available to both
compose and solve matrix equations,
including (a) three options (under- , exact-
and over-relaxation) to estimate the matrix
in the strong coupling loop, (b) three options
(under-, exact- and over-relaxation) to
estimate the matrix in the nonlinear loops
(i.e., the rainfall-evaporation loop in the
subsurface module and the solute transport-
chemical reaction coupling loop in the
reactive chemical transport module), (c)
four options (successive subregion block
iteration, successive point iteration,
polynomial preconditioned conjugate
gradient, and incomplete Cholesky
preconditioned conjugate gradient
methods) to solve the linear/linearized
matrix equations, (d) two options (consistent
and lumping) to treat the mass matrix, (e)
desired time weight factor, ranging from 0.0
to 1.0, to change the type of numerical
scheme, and (f) two options, (the N+1
upstream or the Galerkin weighting
functions) to compose the matrix equations
when transport is simulated by the
conventional Eulerian finite element
method; (13) The time step size can be
reset automatically when boundary
conditions and/or sources/sinks change
abruptly; (14) The examination of mass
balance over the entire domain of interest
for every time step is provided; (15) Many
types of chemical reactions are included in
3DHYDROGEOCHEM, including aqueous
complexation, adsorption/desorption, ion-
exchange, precipitation/dissolution, redox,
and acid-base reactions; (16) The extension
of chemical reactions from chemical
equilibrium to kinetics can be implemented
without much difficulty informulation. Table
1 briefly describes 3DHYDROGEOCHEM
from another view.
DESCRIPTION OF THE MODEL
Mathematical Statements
3DHYDROGEOCHEM is designed to
solve the following system of equations
along with initial and boundary conditions,
which describes water flow, heat transfer,
and reactive chemical transport through
saturated-unsaturated porous media.
The Governing Equation of
Subsurface Flow
The governing equation for flow is
basically the Richards'equation with density
effect taken into account. Under the
assumption of rigid media, the governing
equation can be written as follows.
p d9 3h
-z --- = V
Po
3t
K Vh + —Vz
Po
(1)
where h is the pressure head; t is time; V is
the divergence operator; K is the hydraulic
conductivity tensor; z is the potential head;
q is the source and/or sink; p is the
groundwater density with dissolved
chemical concentrations; po is the
referenced groundwater density; p* is the
density of eitherthe injected orthe withdrawn
groundwater; and 9 is the moisture content.
The hydraulic conductivity K is given by
Table 1. The Property of 3DHYDROGEOCHEM
Model Number
Model
Authors
Year
Dimensions
Coupling Relation a
Coupling Method b
Activity Coefficient °
Aqueous Complexation
Sorption/Desorption d
Precipitation/Dissolution
Redox
Temperature B
18
3DHYDROGEOCHEM
Cheng and Yeh
1995
3
F-H-C-T
Direct and 2-step
Davies
yes
IE, SC, SL, and TL
yes
yes
V
a F-H-C-T denotes coupling among the flow equation, the heat transfer equation, the chemical
relations, and the transport equations.
b Direct denotes strong coupling among the flow equation, the heat transfer equation, the
chemical relations, and the transport equations; 2-step, weak coupling among flow, transfer,
and transport but still strong coupling between the chemical relations and the transport
equations.
c D-H denotes the extended Debye-Huckel equation.
dcIE denotes ion-exchange; SC, surface complexation; SL, single-layer model; and TL, triple-
layer model.
o The temperature distribution can be either read in or computed.
-------
pg
K = —k =
-k k =
-K k
(2)
where |i is the dynamic viscosity of the
groundwater with dissolved chemical
concentrations; po is the referenced
groundwater dynamic viscosity; g is the
gravity; k is the permeability tensor; ks is
the saturated permeability tensor; k. is the
relative permeability or relative hydraulic
conductivity; and Kso is the referenced
saturated hydraulicconductivitytensor. The
referenced values are usually taken at zero
dissolved chemical concentration. The
Darcy flux is calculated as follows.
V = - K •
V h + Vz
(3)
To solve the governing flow equations in
transient-state simulations, the pressure
head over the domain of interest must be
prescribed as the initial condition. This
initial distribution of pressure head can be
obtained by either field measurements or
by solving the steady-state version of
Equation (1). Four types of boundary
conditions, from a mathematical point of
view, are considered to deal with a variety
of physical phenomena that can occur on
the boundary. They are (1) prescribed total
head (Dirichlet) boundary conditions which
are used when the hydraulic head is known
as a function of time on the boundary, (2)
Neumann boundary conditions which are
used when the groundwater flux due to the
pressure head gradient is described as a
function of time on the boundary, (3) Cauchy
boundary conditions which are employed
when the total groundwater flux, due to
both the pressure and the potential head
gradients, is described as afunction of time
on the boundary, and (4) Variable boundary
conditions which are usually applied to the
soil-air interface boundary (see the original
report for full description).
The Governing Equations of
Reactive Chemical Transport
Before the governing equations are
described, thefollowing nomenclature used
in the model needs to be introduced:
(4)
(5)
i=ndps+1
(6)
(7)
(8)
(9)
where T. and T are the total analytical and
modified total analytical concentrations of
the j-th component, respectively; C., S., P.,
and P. are the total dissolved, total sorbed,
total precipitated, and modified total
precipitated concentrations of the j-th
component, respectively; a is the
concentration of the j-th component; x: and
ayx are the i-th complexed species and its
associated stoichiometric coefficient with
respect to the j-th component, respectively;
y. and ayy are the i-th adsorbed species and
its associated stoichiometric coefficient with
respect to the j-th component, respectively;
z. and ayz are the i-th ion-exchanged species
and its associated stoichiometric coefficient
with respect to the j-th component,
respectively; p: and arp are the i-th
precipitated species and its associated
stoichiometric coefficient with respect to
the j-th component, respectively; Mx, My,
Mz, and Mp are the numbers of the
complexed, adsorbed, ion-exchanged, and
precipitated species, respectively; ndps is
the number of the dominating precipitated
species which are precipitated species with
a concentration much higher than the total
dissolved concentrations of their associated
aqueous components. The aqueous
components are mobile in subsurface
systems, whereas the adsorbent
components are immobile. The aqueous
components may be involved in the
complexation, adsorption/desorption, ion-
exchange, and precipitation/dissolution
reactions, whereas the adsorbent
components can only be involved in the
adsorption/desorption reaction. From
another viewpoint, the aqueous
components exist at least in the aqueous
phase, whereasthe adsorbent components
exist only in the solid phase. The complexed,
adsorbed, ion-exchanged, and precipitated
species are the products of the
complexation, adsorption, ion-exchange,
and precipitation reactions, respectively.
The dominating precipitated species are
arranged as the first ndps precipitated
species so that they can be easily indexed
in both derivation and formulation. To
compute reactive chemical transport, three
types of governing equations are basically
solved: the governing equations of aqueous
components, of adsorbent components, and
of ion-exchange sites. Two approaches,
defined as Approach 1 and Approach 2 in
this report, are considered to deal with the
reactive chemical transport of aqueous
components.
For aqueous components:
< Approach 1 >
3C d(s +
9 L+V-VC +9——
dt ' dt
( 36 p' V
+K+—+—* V,
\ dt p )
-v-(eo-vF)-—Fv — •
P ' \P,J
= -v •[QD-V(SI + F)]
,96 P* / -N
} + - Ct + — q(Si + P)
dt ' p
L dt
: Approach 2
(10)
ar
—-+v•^«
d t '
V -(QD-A.,^
1 p 9'
NON
T V -(BD-
I VEC
9i
(11)
-------
3 C
3 T
3 C
3 T
3 C
3EC
(12)
For adsorbent components:
For ion-exchange sites:
je[Na+1,NON]
(13)
)j ie[l,NSITE]
U L U L
(14)
where 6 is the water content; if is the first-
order decay rate of the i-th precipitated
species; EC. is the ion-exchange capacity
of the i-th ion-exchange site; t is time; V is
theDarcyflux; Visthe del operator; Disthe
dispersion coefficient tensor; Kj and >UEC are
the decay constants of the j-th component
and the i-th ion-exchange site, respectively;
3( ) /d t is the partial derivative with
respect to time; C* is the total dissolved
concentration of the j-th component in the
source; Na is the number of aqueous
components; NON and NSITE are the
numbers of components and ion-exchange
sites, respectively; and A^, A^, and AjNON+i
are the concentration' derivatives of
(1
-------
je[Na+1,Na+Ns]
(18)
Mole balance equations for ion-exchange
sites:
NOMZJ(i)+NOMZI(i)
k-NOMZJ(i)+1
ie[l,NSITE]
(19)
Mass action equation for aqueous
complexation:
ie[l,Mx]
(20)
Mass action equation for adsorption/
desorption:
ie[l,My]
(21)
Mass action equation for ion-exchange
reactions:
K
B,.
(22)
Mass action equation for precipitation/
dissolution:
(23)
where Ns is the number of adsorbent
components; ukis the valence of the k-th
ion-exchanged species; NOMZI(i) is the
number of ion-exchanged species involved
in the i-th ion-exchange site; NOMZJ(i) is
the number of ion-exchanged species
involved in the 1 -st through the (i-1 )-th ion-
exchange sites; Aj* is the activity of the i-th
complexed species; K* is the equilibrium
constant of the i-th complexed species; X( is
the activity of the j-th component species;
B* is the activity of the i-th adsorbed species
of surface complexation; K,y is the
equilibrium constant of the i-th adsorbed
species; LNI(i) is the number of the ion-
exchanged species, which is taken as the
reference species for the i-th ion-exchange
site, on the ion-exchanged species list;
K
equilibrium constant of the i-th precipitated
species.
When the dominating precipitated
species exists in the system, we introduce
the modified total analytical concentration,
ratherthan the total analytical concentration,
of components to be the primary dependent
variable forthe system of reactive chemical
transport.
Numerical Strategies
The numerical strategies employed in
3DHYDROGEOCHEM include (1) strong
coupling among flow, heat transfer, and
reactive chemical transport modules for
steady-state simulations to ensure accurate
results, (2) weak coupling amongflow, heat
transfer, and reactive chemical transport
modules for transient-state simulations to
save computer time, (3) Galerkin finite
element method to solve the governing
equation of subsurface flow, (4) (n+1)
upstream weighting finite element to solve
the steady-state version and the
Lagrangian-Eulerian finite element method
to solve the transient-state version of the
governing equations of both reactive
chemical transport and heat transfer, and
(5) Newton-Raphson method to solve the
governing equations of chemical
equilibrium. In solving the governing
equations of reactive chemical transport,
the two approaches employed in the model
to solve for the transport of aqueous
components are emphasized. The
associated numerical strategies, under the
considerations of peak/valley capturing as
well as the existence of dominating
precipitated species, are also described in
the computation dealing with the transport
of aqueous components. The details of all
numerical strategies can be found in
Chapter 3 of the original report.
Sample Problems
The 3DHYDROGEOCHEM model has
been verified with twenty-five example
problems, as stated in the original report. It
also has been applied to solve four three-
dimensional problems, including a
subsurface flow problem, a heat transfer
problem, a reactive chemical transport
problem, and a coupled flow-transport-
transfer problem. The coupled problem is
described in the following.
Problem Description
This sample problem is used to
demonstrate the coupling among
subsurface flow, heattransfer, and reactive
chemical transport. The required universal
parameters for the associated simulation
are listed in Table 2.
For computing the temperature-
dependent groundwater density, the
following equation is employed in this
example problem.
+ a3T2+a4T3
pw(T) = a1+
(24)
where pw(T) is the temperature-dependent
pure water density associated with
temperature T; av a2, a?, and a4 are
temperature-density relationship para-
meters. One thing that should be noted is
that the referenced groundwater density
po(= 997.04882 g/dm3) is associated with
the referenced temperature To (= 298.3 °K).
Two materials are included: material 1 for
the upper half and material 2 for the lower
half as pictured by Figure 1. Material
parameters needed for the simulation are
given in Table 3. In computing subsurface
flow, the following three equations are
employed to describe the soil properties.
K
K
Ksat if h > 0
= K
h, -h
(25)
if h < 0
(26)
Table 2. Universal Parameters Needed for the Simulation of the Sample Problem
is equilibrium constant of the k-th
ion-exchanged species; K.p is the
k.LNI(i)
Universal parameter
Referenced groundwater density, po [g dm3]
Referenced groundwater dynamic viscosity,no [g dm1 day1]
Referenced temperature, To [°K]
Specific heat of groundwater, Cs [dm2 day1 °K1]
Density-temperature relationship parameter 1, a, [dimension/ess]
Density-temperature relationship parameter 2, a2 [dimension/ess]
Density-temperature relationship parameter 3, a3 [dimensionless]
Density-temperature relationship parameter 4, a4 [dimensionless]
Ideal gas constant, R [g dm2 day1 °K1 mole1]
Value
997.04882
7.9056x10s
298.3
3.12035x1015
-410.174
13.2296
-0.0403921
3.97476x105
6.2067x1O15
-------
de
dt
(27)
The temperature-dependent equilibrium
constant is calculated according to the
van't Hoff equation. Table 4 lists chemistry
information needed to compute chemical
equilibriumforasimplif led closed carbonate
system that is considered in this application.
The ionic-strength effect is taken into
account by utilizing the extended Debye-
Huckel formula.
Given well-designed pre-initial and
boundary conditions, the steady-state
simulation can provide a computational
result that presents a uniformly-polluted
domain under an isothermal and nearly
hydrostatic condition. With this result as
the initial condition, the temperature effect
on a removal process can be investigated
by activating a pumping well in the middle
of the domain and controlling the flow-in
groundwater with different temperatures.
The detailed settings are described as
follows. The pre-initial conditions are given
200 dm for the total head, 298.3 °K for the
temperature, and 4x1 Q-4 M for the total
analytical concentrations of both Ca2+ and
CO32~ throughout the entire domain. The
pH-value is kept at 10. The top, bottom,
front, and back boundaries are assumed
impermeable for groundwater, heat, and
chemicals to pass through. To perform the
steady-state solution as mentioned above,
a Dirichlet boundary condition of 200 dm is
specified for the referenced total head on
both the right and left boundaries while
Variable boundary conditions are applied
to those two boundaries for heat transfer
and reactive chemical transport. The flow-
in temperature is 298.3 °K. The flow-in
concentration is 1.2247061 x1Q-4M for both
Ca2+ and CO32-, where 1.2247061x10'4 M is
the total dissolved concentration for the two
components at equilibrium when 4x1 Q-4 M
is the total analytical component
concentration. Due to the density effect, a
truly hydrostatic situation cannot be
achieved. However, a nearly hydrostatic
condition can be performed because the
density effect is not important.
When the transient-state simulation
begins, three point sinks are used to
simulate a pumping well for a removal
purpose. These three point sinks, with a
pumping rate of 5x103 drrWday for each,
are at (250, 100, 80), (250, 100, 90), and
(250, 100, 95). It is thus imaginable that
groundwater will flow into the domain
through the left and the right boundaries.
Relative clean groundwater, with a total
dissolved concentration of 1Q-8 M for
components 1 and 2, is assumed to supply
Figure 1. The simulated domain of the
sample problem,
this flow-in groundwater. A temperature of
308.3 °K is associated with the flow-in
groundwater through the left boundary,
whereas the flow-in groundwater through
the right boundary is still under the
temperature of 298.3 °K. To put a focus on
the temperature effect, the referenced total
head of 200 dm is kept as in the steady-
state simulation, so that the subsurface
flow will be almost symmetric with respect
to x = 250 dm (i.e., the middle of the domain
in the x-direction). It is thus adequate,
directly based on the computational results,
to examine how significantly temperature
can affect reactive chemical transport.
Forty time steps are included in the
transient-state simulation, where 1, 2, and
4 days are set to be time-step sizes of the
first, second, and third through fortieth time
steps, respectively. 5x1 Q-5 dm is assigned
as the absolute error tolerance for a
convergent solution of subsurface flow. In
the meantime, 10~3 and 5x10~5 are used to
be the relative error tolerances in computing
for chemical concentrations and
temperature, respectively.
Simulation Results
As mentioned above, the steady-state
solution is controlled to be almost identical
to the pre-initial condition. That is, the
simulated domain has chemical
concentrations uniformly distributed under
an isothermal and a nearly hydrostatic
condition at the beginning of the transient-
state simulation.
With a pumping rate of 5x103 drrrVday
assigned to the three point sinks in the
middle of the domain, aflowfield, symmetric
in the x-direction with respect to x= 100dm,
Table 3. Material Parameters Needed for the Simulation in the Sample Problem
Material parameter Material 1
Saturated K^ [dm day1]
Saturated K [dm day1]
Saturated Kzz [dm day1]
Saturated Kx (or KJ [dm day1]
Saturated KX2 (or KJ [dm day1]
Saturated Kyz (or Kyz) [dm day1]
Effective porosity, 9gff [dimension/ess]
Soil property parameter 1, 9a [dimension/ess]
Soil property parameter 2, ha [dm]
Specific heat of dry medium, Cm [dm2 day1 °K1]
Apparent thermal conductivity, ~k [g dm day3 °K1]
Bulk density, pb [g dm3]
Longitudinal dispersivity, aL [dm]
Transverse dispersivity, aT [dm]
Diffusion coefficient, am [dm2 day1]
Tortuosity, i [dimension/ess]
1.0
1.0
1.0
0.0
0.0
0.0
0.3
0.15
-1000.0
1.0x10"
2.0x1 013
1200.0
1.0
0.1
0.0001
1.0
Material 2
5.0
5.0
5.0
0.0
0.0
0.0
0.2
0.1
-1000.0
2.0x10"
2.0x1 013
1200.0
2.0
0.2
0.0002
1.0
Table 4. Chemistry Data for the Simulation in Application 4
Product species
Stoichiometric coefficient Log K
Ca2+ CO,2 H+
Reaction enthalpy*
OH
HCO3
H2CO3
Ca(OH)2(s)
CaC°3(s)
0
0
0
1
1
0
1
1
0
1
-1
1
2
-2
0
-14.000
10.200
16.500
-21.900
8.300
4.1 67690x1 O19
-1.112280x1019
-1.1 71 028x1 020
1.256353x1019
9. 353595x1 'O13
* The unit of reaction enthalpy is [g dm2 day2 mole1], 10 [kJ mole1] = 7.46496x1 Ow [g dm2 day2
mole1]
-------
can be observed when identical boundary
conditions are applied to both the left and
the right boundaries. In this application,
however, a temperature of 308.3 °K is given
to the incoming subsurface flow from the
left boundary, which is 10 °K higher than
that initially set for the simulated domain or
that of the incoming subsurface flow from
the right boundary. Thus, the groundwater
density will not be symmetrically distributed
inthex-direction. Based on the fact that the
possible highest temperature is 308.3 °K
and the possible maximum total analytical
concentrations are 4x10'4 M for components
1 (i.e., Ca2+), 2 (i.e., CO/-) and less than
10'4M for component 3 (i.e., H+) through the
whole simulation, the influence of
temperature and chemical concentration
on the groundwater density can be as large
as 0.35% (according to Figure 4.24) and
0.004% (according to Eq. 2.35(b)) during
the transient-state simulation.
Figures 2 and 3 present the total head
distributions at time = 75 days and 155 days,
respectively. The (a) part of these two
figures puts a focus on planes y = 100 dm,
y = 200 dm, and z = 100 dm, while part
(b) concentrates on the planes of z = 50
dm,z=100dm,z=150dm,andy= 100dm.
A symmetry in the y-direction and an
asymmetry in the z-direction can be seen
as expected. Although the symmetry in the
x-direction looks true, it does not exist
actually due to the groundwater density
effect mainly caused by a higher incoming
temperaturefrom the left boundary. Figure 4
plots the flow rate on both the left and the
right boundaries at time = 75 days (part (a))
and 155 days (part (b)). Again, the
magnitudes of the flow rates on the two
boundaries cannot be distinguished
visually. However, the computational result
does provide values to tell the difference.
On the other hand, since the groundwater
density effect is so minor in this case that
the difference mentioned above might be
ignored, the subsurface flow could be
assumed symmetric in the x-direction. The
temperature-dependent removal process
will be discussed later under this
assumption.
Figure 5 shows the temperature
distribution at time = 75 days with emphases
on (a) plane z = 100 dm and (b) plane y =
100 dm. Likewise, Figure 6 plots the
temperature distribution at time = 155 days.
From Figure 5, only the isothermal line of
299 °K implicitly indicates a possible
existence of pumping wells. The pumping,
nevertheless, becomes obvious in Figure 6
because a longer simulation time allows
heat to be transferred from the left boundary
to the neighborhood of pumping wells,
where the flowfield is greatly dominated by
pumping.
Figures 7 and 8 depict the computational
results of simulating a removal process by
pumping at time = 75 days and 155 days,
respectively. In these two figures, part (a)
is used to emphasize the total analytical
concentration distribution on the planes of
z = 65 dm and 135 dm, while part (b) is to
describe the distribution of the total
analytical concentration on planes y = 0 dm
and y = 100 dm. Due to the fact that
Ca(OH) (s) does not exist under the
competition with CaCO3(s) and dissolution is
the only reaction of mass transfer between
the solid and the aqueous phases through
the transient-state simulation, the total
analytical concentrations of component Ca2+
(b)
(b)
Figure 2. The total head distribution at
time = 75 days.
Figure 3. The total head distribution at
time =155 days .
Figure 4. The total head distribution and
boundary flow rates at (a) time
= 75 days and (b) time =155
days.
-------
(b)
Figure 5. The temperature distribution at
time = 75 days.
Figure 6. The temperature distribution at Figure 7.
time =155 days.
The total analytical concentra-
tion distribution of component
1 (Ca2+) or2 (CO/) at time = 75
days.
and CO32- are equal. This can be examined
based on the prescribed equilibrium
constants for such a simplified system.
Therefore, the concentration distributions
shown in Figures 7 and 8 can be for
component Ca2+ or CO32-.
Both figures point out that the higher
incomingtemperaturefromthe left boundary
drives the 0.0003957 M concentration line
on the left side to move faster than that on
the right side. But the 4.65x10'5 M
concentration lines on both sides seem to
move with the same progress toward the
locations of sinks. As a matter of fact, the
higher temperature makes this 4.65x1 Q-5 M
concentration line on the left move slower
than that on the right. Figure 9 provides a
two-dimensional concentration contour plot
on the plane of y = 100 dm, which gives a
clearer view on the two 4.65x1 Q-5 M
concentration lines at time = 155 days. The
contrast between the 0.0003957 M and the
4.65x1 Q-5 M concentration lines can be
explained as follows.
According to the van't Hoff equation, a
higher temperature will increase the
equilibrium constant of a chemical reaction
if the associated reaction enthalpy is
positive. In this application, thus, a higher
temperature coming in from the left
boundary will help the precipitation of
CaCO3(s) and increase the effort of
completely removing the precipitant. This
is why the 4.65x10~5 M concentration line on
the left side migrates slower than that on
the right side.
Given the total analytical component
concentration, on the other hand, a higher
temperature will cause a lower dissolved
component concentration to be achieved at
equilibrium. This lower dissolved
component concentration at an upstream
location andthecurrenttime will bebrought
to a downstream location along with the
subsurface flow at the nexttime. Therefore,
the total analytical component concentration
at that downstream location will be reduced
at the next time step if the temperature of
the upstream location is higher than that of
the downstream location at the current time.
This basically explains what happens to the
0.0003957 M concentration line in this case.
The faster moving concentration line on the
left part of the domain causes less dissolved
calcium (or carbonate) to migrate from
upstream to downstream due to a higher
upstream temperature.
The above discussion, in away, implicitly
signifies the need of numerical models for
quantitative analyses: a highertemperature
introduced from the upstream boundary will
cause more precipitation of CaCO3(s) within
the upstream region due to a positive
reaction enthalpy, but it also generates a
continuous dissolution within the
downstream region due to the reduced total
analytical concentration (as explained
above) at a lower temperature. Under a
lower temperature, a downstream location
always has less total dissolved
concentration coming in than going out,
which makes this location undersaturated
with the existence of precipitated species.
To attain equilibrium, dissolution occurs to
decrease the concentration of CaCO3(s) at
this downstream location. These opposite
situations hence lead to a difficulty in
determining whether or not an introduction
of a higher temperature from the upstream
boundary would reduce the overall removal
effort. To this point, a qualified numerical
-------
model is required. One should be aware
that the upstream and the downstream
regions mentioned above are actually
representing a relatively high-temperature
and a relatively low-temperature region.
Therefore, their sizes are to be changed
with the transfer of heat. The speed of heat
transfer, relative to that of reactive chemical
transport, has been found able to
significantly affect the removal result. This
subject, even though not discussed here,
could be interesting to environmental
researchers and engineers.
Figure 8. The total analytical concentra-
tion distribution of component
1 (Ca2+) or 2 (CO/-) at time
=155 days.
200
175
150
125
100
75
50
25
0.0003957
50
100 150 200 250 300 350 400 450
500
Figure 9. The total analytical concentration distribution on the plane of y = 100 dm for
component 1 or 2 at time = 155 days.
-------
Gour-Tsyh (George) Yeh and Hwai-Ping (Pearce) Cheng are with the Pennsylvania
State University, University Park, CA 16802.
Thomas £. Short is the EPA Project Officer (see below).
The complete report, entitled "3DHYDROGEOCHEM: A 3-DimensionalModel of Density-
Dependent Subsurface Flow and Thermal Multispecies-Multicomponent
HYDROGEOCHEMical Transport, "(OrderNo. PB99-150419; Cost:$41.00, subject
to change) will be available only from:
National Technical Information Service
5285 Port Royal Road
Springfield, VA22161
Telephone: 703-487-4650
The EPA Project Officer can be contacted at:
U.S. Environmental Protection Agency
National Risk Management Research Laboratory
Subsurface Protection and Remediation Division
P.O. Box 1198
Ada, OK74820
United States
Environmental Protection Agency
Center for Environmental Research Information
Cincinnati, OH 45268
BULK RATE
POSTAGE & FEES PAID
EPA
PERMIT No. G-35
Official Business
Penalty for Private Use
$300
EPA/600/SR-98/159
------- |