Simulation of bimolecular and trimolecular reactive transport in porous ...
WATER RESOURCES RESEARCH, VOL. ???, XXXX, DOI:10.1029/,
1
2
3
Simulation of bimolecular and trimolecular reactive transport in porous media using colocation probability functions. 1
1
J. M. Barnard, C. E. Augarde,
C. E. Augarde, School of Engineering and Computing Sciences, Durham University, South Road, Durham, DH1 3LE, UK. ([email protected]) J. M. Barnard, School of Engineering and Computing Sciences, Durham University, South Road, Durham, DH1 3LE, UK. ([email protected]) 1
School of Engineering and Computing
Sciences, Durham University, South Road, Durham, DH1 3LE, UK.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 2 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 4
Abstract.
Random walk particle tracking methods have shown them-
5
selves to be effective for describing mass transport in porous media due to
6
their ability to completely remove the problem of numerical dispersion. Sim-
7
ulation of reactions between different chemical species is a much less advanced
8
topic than in other methods of simulating mass transport in porous media
9
and has so far been restricted to simple reactions of the form A + B → C.
10
Here a method for simulating reactions with a stochastic particle tracking
11
framework is developed and discussed with respect to previous methods. This
12
method is based on colocation probability functions and is capable of deter-
13
mining the likelihood of the colocation of any number of particles with dif-
14
ferent arbitrary, Gaussian location probability functions. The method is val-
15
idated against experimental data from Gramling et al. [2002] and then used
16
in an expanded A + B + C → P form as a proof of concept for larger, more
17
complicated systems.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X-3
1. Introduction 18
The accurate simulation of reactive mass transport through porous media is an im-
19
portant capability for the study of various different environmental phenomena such as
20
contaminant transport within soil. As part of the multi-disciplinary, ROBUST project
21
(Regeneration Of Brownfield Using Sustainable Technology), which aims to develop a new
22
remediation technologies for contaminated land a new computational model for reactive
23
mass transport in porous media has been developed using a particle tracking based method
24
for simulating transport in tandem with a finite element porous media flow model.
25
26
The physical basis for porous media transport models is generally the advectiondispersion-reaction equation (ADRE) [Bear and Cheng, 2010];
27
∂C + v · ∇C = ∇ · (D∇C) + Q ∂t
28
where C is solute concentration, t is time, v is pore water velocity, D is a combined
29
dispersion-diffusion parameter, Q represents all of the sources and sinks for chemical
30
species including reactions and ∇ is the differential operator. There are many studies
31
into reactive mass transport through porous media that use either analytical or numerical
32
solutions to the ADRE as the basis of their work [Berkowitz et al., 2006; Edery et al.,
33
2010; Gramling et al., 2002].
(1)
34
Most models used to simulate porous media flow use a continuum approach to describe
35
the medium and its pore fluids. This works well for fluid flow problems but once solute
36
transport is added to the problem it becomes unviable for a range of reasons. Describing
37
the transport of a solute concentration field using a continuum based approach such as the
38
finite element method leads to smearing or numerical dispersion of the concentration field
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 4 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 39
parallel to the mean flow vector particularly in areas with high concentration gradients.
40
Though the extent of this effect can be reduced by refining the spatial resolution of the
41
simulation this is an inefficient and inelegant way to deal with the problem.
42
Another disadvantage of the continuum approach to describing mass transport in porous
43
media is that it cannot describe small scale variation in the solute concentration field. It
44
therefore assumes that reactants are mixed perfectly on the element scale. Experiment
45
[Gramling et al., 2002] and pore scale simulation [Cao and Kitanidis, 1998] have both
46
shown that the assumption of perfect mixing overpredicts the amount of a reactant pro-
47
duced.
48
Particle tracking techniques such as the discrete time random walk (DTRW) and con-
49
tinuous time random walk (CTRW) techniques completely eliminate the problem of nu-
50
merical dispersion from simulations of the ADRE and have been used in a wide range
51
of fields [Montroll and Weiss, 1965; Prickett et al., 1981]. Other advantages of particle
52
tracking methods are that they have been shown to be adept at describing anomalous
53
or non-Fickian transport [Berkowitz et al., 2006] and imperfect mixing of reactants in a
54
porous medium [Edery et al., 2010].
55
The application of the DTRW to the solution of the ADRE as set out by Prickett et al.
56
[1981] describes the concentration field of a chemical species using a field of particles.
57
These particles behave independently of all other particles of the same species and can be
58
located at any point within the domain.
59
In the DTRW method particles are advected according to the porous media flow field
60
and then walked a random distance ∆L over a fixed time step ∆t. The particles can take
61
any location within the domain unlike simulations using the CTRW method in which
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X-5
62
particle locations are limited to a lattice of locations. Simulating reactive mass transport
63
using either of the random walk methods is a more convoluted process than in any of the
64
continuum based approaches. As particles are treated as indivisable units of solute they
65
are reacted on a particle to particle basis. For example in the simulation of a reaction of
66
the form A + B → C, if a reaction takes place then an A and a B particle are destroyed
67
and a new C particle is formed.
2. Methods 68
In the random walk method each particle represents a ‘packet’ of the chemical species
69
rather than an individual particle, but is treated as though this is the case. The packet of a
70
chemical species discribed by a particle cannot change shape due to the effects of dispersion
71
or react partially with a packet of another species. These particles are advected in the
72
direction of the flow field and then dispersed a random distance parallel and perpendicular
73
to the direction of flow. Particle reactions are then performed on a particle by particle
74
basis.
75
The evolution of the particle field is performed in three stages; advection, dispersion
76
and reaction. The advection stage is performed by moving each particle according the
77
average flow field through the porous medium. This average flow field is calculated from
78
nodal pressure values using the methods of Cordes and Kinzelbach [1992] and Pollock
79
[1988]. These nodal pressure values could be determined using any continuum based
80
approach describing porous media flow. Here the finite element method is used with
81
bilinear quadrilateral elements.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 6 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
2.1. Discrete Time Random Walk Method 82
In the Pollock method particles are moved according to the flow field at their initial
83
location. The explanation of the method, from Pollock [1988], is for streamlines in two
84
dimensions. The expansion to three dimensions is trivial and so will not be given here.
85
The x and y components of a particle’s velocity are given as;
86
vxp =
(2)
vyp
(3)
87
1 (vx − vx1 )(xp − x1 ) + vx1 ∆x 2 1 = (vy − vy1 )(yp − y1 ) + vy1 ∆y 2
88
where vxp is the x component of the velocity of the flow field at the particle location,
89
∆x is the width of the element in which the particle is currently located, vx1 and vx2
90
are the x components of the flow field velocity at x1 and x2 , respectively, of the element,
91
xp is the x coordinate of the particle location and x1 and x2 are x coordinates of the
92
element boundaries perpendicular to x. The y components and coordinates of the above
93
parameters are labelled accordingly. This method was initially developed for use with
94
the finite difference method. It can be used unaltered with the finite element method if
95
regular bilinear quadrilateral elements are being used. If isoparametric elements are being
96
used then the global coordinates must be transformed so the element is square in the local
97
coordinate system.
98
The differential with respect to time of the x component of particle velocity is; dvx dt
99
!
= p
dvx dx
!
dx dt
!
.
(4)
p
100
The differential of vx with respect to x = (vx2 −vx1 )/∆x and the change of particle position
101
with respect to time is the definition of vxp . Therefore the simplification of (4) gives; dvx dt
102
D R A F T
!
= p
1 (vx − vx1 )vxp ∆x 2
May 31, 2012, 12:44pm
(5)
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 103
which can be rearranged to give; 1 1 dvxp = (vx − vx1 )dt. vxp ∆x 2
104
105
X-7
(6)
Integrating (6) between times t1 and t2 (t2 > t1 ) gives, "
#
106
vxp (t2 ) 1 ln = (vx − vx1 )∆t vxp (t2 ) ∆x 2
107
where ∆t = t2 − t1 . Rearrangement of (7) gives travel time from any x position within an
108
element to any other as; ∆tx =
109
110
∆x vx2 − vx1
!
"
(7)
#
vx2 ln . vx1
(8)
Analogous statements can be made for a particle’s movement with respect to y.
111
The first step in calculating the movement of a particle is to determine which element
112
edge it will leave through. This is done by calculating the travel times between a particle’s
113
location and each element edge using (8). Assuming the particle is not on an element
114
boundary the smallest positive travel time will correspond to the edge through which the
115
particle will exit. If the particle is on an element boundary then the travel time to that
116
boundary is obviously zero but the particle should only exit through that element edge if
117
the flow direction is towards it. To address this, if a particle is on an element boundary,
118
the element which it is in is determined by the direction of the pore fluid flux at that
119
boundary, the particle being placed in the element it is about to advect into rather than
120
advect out of. In addition if a particle is on an element edge, travel times of zero are
121
ignored when determining which edge the particle will exit through.
122
Dispersive transport is then simulated by the random walk of the particles. Each particle
123
is moved a random distance in the directions parallel and perpendicular to the mean flow
124
vector along which the particle was advected. A pseudo random number generator with D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 8 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 125
an N[0,1] distribution produced using a Box-M¨ uller transform on a U[-1,1] pseudo random
126
number generator is used. The generator is then scaled to the dispersivity and the time
127
step. The location probability function (LPF) for each particle is therefore Gaussian in
128
form with a standard deviation, σ, dependent on dispersivity and time step length. The
129
time dependence of dispersive transport can take one of two forms, either the Fickian form
130
where σ = (Dh ∆t)−0.5 or a truncated power law form (TPL) [Edery et al., 2010; Berkowitz
131
et al., 2006] where σ = e−t/t2 (1 + t/t1 )−1−β where β is the power law exponent, t1 is the
132
characteristic transition time and t2 is the cutoff time to Fickian transport. This TPL
133
form of time dependency was used with the CTRW method by Dentz et al. [2004], hence t1
134
being labelled in accordance with CTRW terminology. The TPL form of time dependency
135
has been used due to its ability to describe anomalous or non-Fickian transport which
136
occurs in heterogenous porous media [Berkowitz et al., 2000; Dentz et al., 2004; Edery
137
et al., 2010]. Edery et al. [2010] have shown that use of the TPL form of time dependency
138
can produce closer fits to experimental data for random walk based simulations.
139
Reactions in the random walk methods can be handled in several ways. Edery et al.
140
[2010] react particles if their separation is less than a reaction radius R. This reaction
141
radius is set empirically to fit the experimental data. Benson and Meerschaert [2008]
142
determine the likelihood that particles will be colocated at the end of the time step as a
143
function of current separation and time step size, or that for the next time step relative
144
motion will equal current separation. If this probability is greater than the value of a
145
random number from a U[0,1] generator then the particles are reacted. The probability of
146
this occurance is given by the convolution of two Gaussian LPFs which is also a Gaussian
147
function. Benson and Meerschaert [2008] approximate this Gaussian colocation probabil-
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X-9
148
ity function (CLPF) with a tent function to increase computational efficiency. This works
149
well when the location probability functions only have a single value for σ. If the LPFs
150
are of the more complicated form, such as when longitudinal and transverse dispersivity
151
are not equal, then determining the properties of the CLPF and the tent function to be
152
used become more complicated. This reduces any computational advantage gained from
153
the use of a tent function. 2.2. Particle Based Reaction Methods
154
In this study reaction is handled in a stochastic fashion based on the probability of one
155
particle of each of the reactant species being colocated over a time step ∆t and is similar
156
in some respects to the method used by Benson and Meerschaert [2008]. A particle’s
157
probability location function is of the same form as the distribution of the pseudo random
158
number generator used in the dispersive transport step. The probability of a pair of
159
particles meeting at a certain point will be the convolution of each of their LPFs. As the
160
LPFs have Gaussian forms their convolution is trivial. The complete definite integral of
161
this convolution or colocation probability function will then give the probability that the
162
particles will meet at any location during the time step. This method will initially be
163
demonstrated for the convolution of two arbitraty Gaussian functions in two dimensions.
164
Expansion to the convolution of any number of arbitrary Gaussian functions will then be
165
shown to be trivial.
166
167
The location probability function of a particle is expressed as an arbitrary two dimensional Gaussian function with its long axis parallel to the mean flow vector; 2 +2b(x−x
z = He−(a(x−x0 )
168
D R A F T
0 )(y−y0 )+c(y−y0 )
May 31, 2012, 12:44pm
2)
(9)
D R A F T
X - 10 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 169
where; a=
170
cos2 θ sin2 θ + 2σx2 2σy2
(10)
sin2θ cos2θ + 4σx2 4σy2
(11)
sin2 θ cos2 θ + 2σx2 2σy2
(12)
171
b=−
172
173
c=
174
175
and H is the maximum value of z located at (x0 , y0 ) and θ is the angle of the mean flow
176
vector from the y axis..
177
The convolution of two of these functions is given by;
178
z = H1 H2 exp(−(a1 (x − x1 )2 + a2 (x − x2 )2 2b1 (x − x1 )(y − y1 ) +2b2 (x − x2 )(y − y2 ) + c1 (y − y1 )2 + c2 (y − y2 )2 ))
179
180
which can be rearranged to give; 2 +by 2 +cxy+dx+f y+g
z = Heax
181
182
(13)
(14)
where;
183
a = (a1 + a2 ),
(15)
b = (c1 + c2 ),
(16)
c = 2(b1 + b2 ),
(17)
d = −2(a1 x1 + a2 x2 + b1 y1 + b2 y2 ),
(18)
f = −2(c1 y1 + c2 y2 + b1 x1 + b2 x2 ),
(19)
g = a1 x21 + a2 x22 + c1 y12 + c2 y22 + 2b1 x1 y1 + 2b2 x2 y2 ,
(20)
184
185 186
187 188
189 190
191 192
193
194
and
195
D R A F T
H = H1 H2 . May 31, 2012, 12:44pm
(21) D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 11
196
197
The definite integral of the above function is equal to; 1
2Hπ(4ab − c2 )− 2 e
198
−4abg+af 2 +bd2 +c2 g−cdf 4ab−c2
(22)
199
where, for the convolution of n arbitrary Gaussian functions each centred at (xn , yn ), the
200
coefficients a, b, c, d, f and g are given by; a=
201
202
b=
203
204
c=
205
n X
cos2 θi sin2 θi + 2 2σy2i i=1 2σxi
n X
sin2 θi cos2 θi + 2 2σy2i i=1 2σxi
n X i=1
206
d = −2
207
n X
"
i=1 208
f = −2
209
n X
"
i=1 210
211
g=
n X i=1
"
(23)
−
(24)
sin2θi sin2θi + 2σx2i 2σy2i
(25)
cos2 θi sin2 θi sin2θi sin2θi + xi + − + yi 2 2 2σxi 2σyi 4σx2i 4σy2i
(26)
sin2θi sin2θi sin2 θi cos2 θi + y + − + xi i 2σx2i 2σy2i 4σx2i 4σy2i
(27)
#
#
"
#
"
#
cos2 θi sin2 θi 2 sin2 θi cos2 θi 2 sin2θi sin2θi + x + + yi + 2 − + xi yi (28) i 2 2 2 2 2σxi 2σyi 2σxi 2σyi 4σx2i 4σy2i #
"
#
"
#
212
The maximum value of the CLPF, H, is set to one so that the maximum probability of
213
colocation is one and this occurs when the particles are colocated at the beginning of the
214
time step.
3. Results & Discussion 215
The model outlined above will be verified against experimental data from Gramling
216
et al. [2002] showing the transport and reaction of two solutes to form a thrid. The
217
experiment was performed using a box 36 cm long, 5.5 cm high and 1.8 cm deep filled
218
with cryolite sand. The non-reactive transport experiment starts with the sand saturated
219
with a 0.02 M solution of sodium ethylenediaminetetraacetic acid (NaEDTA2− ) which
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 12 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 220
is then displaced by a 0.02 M solution of CuSO4 at a flow rate of 2.67 mL/min, which
221
gives an average pore fluid velocity of 1.21 x 10−2 cm/s. The formation of the product
222
CuEDTA2− was then measured [Figure 1].
223
The experiment is simulated as a two dimensional domain 36 cm long and 5.5 cm high.
224
The flow field exists at a steady state producing a pore fluid velocity of 1.21 x 10−2 cm/s.
225
The left and right hand side boundaries are set as fixed pressure boundaries and top and
226
bottom boundaries are set as no flow boundaries. This initial concentration C 0 of each
227
reactant corresponds to a particle density of 225 per cm2 .
228
The simulation outputs correspond well to the experimental results from Gramling
229
et al. [2002] [Figures 2, 3 and 4]. The shape and peak concentration both match the
230
experimental data well. To produce these results the reaction probability has not been
231
empirically fitted to the experiment like the method used by Edery et al. [2010]. For a
232
Fortran 90 implementation of the model the run time to reach the situation shown in
233
Figure 4 with a time step length of 1 s on notebook computer with a 2.10 GHz, dual core
234
processor and 3 GB of RAM was 1505 seconds. Note that due to the probabilistic nature
235
of the model runtimes are variable. Run times are with 10% of those for simulations using
236
reaction radius method in an otherwise identical code.
237
The overall process of the reaction algorithm is; select particle A, determine which
238
B particles are within a set radius R of the selected A, loop over the set of B and for
239
each calculate (22) and compare to the output a U[0,1] pseudorandom number generator.
240
The checking of particles against one another to determine the probability of a reaction
241
occuring is the most computationally expensive part of the simulation and so certain
242
tactics were implemented within the code to attempt to minimise this. The set of B
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 13
243
less than R away from the selected A are placed in proximity order and so also order of
244
decreasing probability of reaction. In addition if the probability of a reaction occurring
245
ever drops below 10−3 then the loop over the set of B is terminated. Together these
246
additions to the code reduce the number of times that (22) is calculated. The checking of
247
each B particle to see if it is within R of the selected A is computationally expensive if the
248
entire set of B particles is checked. Instead the set of B particles is divided into subsets of
249
the particles within a set distance of each node at the beginning of the subroutine. The
250
closest node to the selected A is then determined and only the subset of B associated with
251
the same node are checked against.
252
Figures 5 & 6 show the variation of the solution at t = 916s with varying time step.
253
For both ∆t = 1s [Figure 3] and ∆t = 2s [Figure 6] the fit to the experiment data is
254
good, though at ∆t = 0.5s [Figure 5] the simulation slightly over-predicts the amount
255
of CuEDTA produced. The principal advantage of the CLPF method is that, for the
256
experiment simulated here, it does not need fitting empirically, unlike the reaction ra-
257
dius method. It scales automatically with the diffusivity and dispersivity of the situation
258
including situations with different logitudinal and transverse dispersivities, thereby au-
259
tomatically adjusting itself to variable flow fields. The stated advantage of both the
260
reaction radius and tent function methods, that they reduce run times by removing the
261
computationally expensive parts of the CLPF method has been shown to be minimal.
262
The second simulation run is used as a proof of concept for reaction involving a higher
263
number of reactants, in this case the reaction is of the form A + B + C → D. The basic
264
setup is the same as Figure 1 with a domain filled with a solution of one of the reactants
265
and the other two reactants flowing in at the same concentration. The setup is altered so
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 14 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 266
that the domain is 20 cm high and 30 cm long and the inflow of reactants into the domain
267
is different. The other reactants then flow into the domain under a constant flow regime
268
from the left hand side. One reactant enters from the top fifth of the left hand side and
269
the third from the lowest fifth. The average pore fluid velocity is 1.21 x 10−1 cm/s and the
270
diffusivities of the reactants and products are all 1.75 x 10−2 . The initial concentration
271
of A in the domain is 900 particles per cm2 Figures 7 & 8 show concentrations of D at
272
100s and 200s respectively. The concentration field is a Gaussian surface weighted by
273
particle locations. Product concentrations relative to reaction concentrations are lower
274
than those in previous example due to the much lower probability of the three particles
275
being colocated than two.
276
Benson and Meerschaert [2008] state that for the sake of computational efficiency re-
277
actions of the type A + B + C → P should be handled as two consecutive reactions, A
278
+ B → D and C + D → P. Whilst it is true that the simulation of an A + B + C → P
279
reaction is computationally demanding it is not necessarily true than a different method
280
should be used. By splitting A + B + C → P into A + B → D and C + D → P the
281
eventuality of D particles existing unreacted in the domain becomes a possibility. If this
282
is possible in the chemical reaction being simulated then it is not a problem, but if the
283
reaction can only take in A,B and C and produce P then it becomes a problem. Also
284
the probability of three particles being colocated is much less than the probability of two
285
particles being colocated, reducing the reaction rate. Which reaction rate is the true one
286
would then affect which would be the better method to use.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 15
4. Conclusion 287
The colocation probability function method for simulating reaction in particle tracking
288
based models of reactive mass transport in porous media flow has been expanded from
289
the version presented by Benson and Meerschaert [2008]. This expansion has included
290
the capability to simulate reactions between any number of particles, each with different,
291
arbitrary Gaussian probability location functions. This method has also been tested
292
against experimental data from Gramling et al. [2002] where it appears to perform well.
293
The simulation does not require the variation of any empirical factors to ensure a good
294
fit to the experiment. As a proof of concept, the multicomponent reaction method has
295
been tested on a reaction of the form A + B + C → D. The overall purpose of this
296
work has been to develop a method for simulating reaction in particle tracking based
297
models of transport in porous media that is capable of simulating any of the reactions
298
that take place within the subsurface. The expansion of this method to include other
299
capabilities such as the simulation of reversible reactions, radioactive solutes, sorbtion of
300
solutes to the soil skeleton, catalysed and non-instantaneous reactions is trivial from a
301
programming perspective. Therefore a framework which is capable of describing a wide
302
range of chemical processes and so a large proportion of the processes involved within soil
303
chemistry has been developed.
304
305
Acknowledgments. This research is part of the Regenerating Brownfield Using Sustainable Technology project. EPSRC Challenging Engineering Grant: EP/G028958/1.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 16 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
References 306
J. Bear and A. H.-D. Cheng. Modeling Groundwater Flow and Contaminant Transport,
307
volume 23 of Theory and Application of Transport in Porous Media. Springer, 2010.
308
D.A. Benson and M.M. Meerschaert. Simulation of chemical reaction via particle track-
309
ing: Diffusion-limited versus thermodynamic rate-limited regimes. Water Resources
310
Research, 44:W12201, 2008.
311
B. Berkowitz, H. Scher, and S.E. Silliman. Anomalous transport in laboratory-scale,
312
heterogeneous porous media. Water Resources Research, 36(1):149–158, 2000. ISSN
313
0043-1397.
314
B. Berkowitz, A. Cortis, M. Dentz, and H. Scher. Modeling non-Fickian transport in
315
geological formations as a continuous time random walk. Reviews of Geophysics, 44(2),
316
2006. ISSN 8755-1209.
317
318
J. Cao and P. K. Kitanidis. Pore-scale dilution of conservative solutes: An example. Water Resources Research, 34:1941–1949, 1998.
319
C. Cordes and W. Kinzelbach. Continuous groundwater velocity fields and path lines in
320
linear, bilinear, and trilinear finite elements. Water resources research, 28(11):2903–
321
2911, 1992. ISSN 0043-1397.
322
M. Dentz, A. Cortis, H. Scher, and B. Berkowitz. Time behavior of solute transport in
323
heterogeneous media: Transition from anomalous to normal transport. Advances in
324
Water Resources, 27(2):155–173, 2004.
325
326
Y. Edery, H. Scher, and B. Berkowitz. Particle tracking model of bimolecular reactive transport in porous media. Water Resources Research, 46(7):W07524, 2010.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 17
327
C.M. Gramling, C.F. Harvey, and L.C. Meigs. Reactive transport in porous media: A
328
comparison of model prediction with laboratory visualization. Environmental science
329
& technology, 36(11):2508–2514, 2002.
330
331
332
333
E.W. Montroll and G.H. Weiss. Random walks on lattices. II. Journal of Mathematical Physics, 6:167, 1965. D.W. Pollock. Semianalytical computation of path lines for finite-difference models. Ground Water, 26(6):743–750, 1988.
334
T.A. Prickett, T.G. Naymik, and C.G. Lonnquist. A’ Random Walk’ Solute Transport
335
Model for Selected Groundwater Quality Evaluations. Illinois State Water Survey Bul-
336
letin, 65, 1981.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 18
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
36 cm
0.02 M
NaEDTA2-
2.67 mL/min
Figure 1. i)
5.5 cm
0.02 M
CuSO4
Schematic diagram showing the experimental setup used by Gramling et al. [2002] t=n
ii)
t = n+1
and then simulated using the DTRW method. L
0.45
C
C
0.4
CuEDTA/C0
0.35
iii)
Simulation Experiment
x
0.3
t=n
iv)
t = n+1
0.25 0.2 0.15 0.1
C
x
L
C
0.05 0 0
x5
10
15
20
x
x (cm)
Figure 2. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 619 s.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 19
0.45 Simulation Experiment
0.4
CuEDTA/C0
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0
5
10
15
20
x (cm)
Figure 3. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 916 s.
0.45 Simulation Experiment
0.4
CuEDTA/C0
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 5
10
15
20
25
x (cm)
Figure 4. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 1510 s.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 20
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
0.45 Simulation Experiment
0.40
CuEDTA/C0
0.35 0.30
0.25 0.20 0.15
0.10 0.05 0.00 0
5
10
15
20
x (cm)
Figure 5. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 916 s for ∆t = 0.5s.
0.45 Series1 Series2
0.4
CuEDTA/C0
0.35 0.3
0.25 0.2 0.15
0.1 0.05 0 0
5
10
15
20
x (cm)
Figure 6. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 916 s for ∆t = 2s.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
Figure 7.
X - 21
Pseudocolour plot of product concentration relative to maximum reactant concen-
tration at t = 100s. Flow direction is from left to right.
Figure 8.
Pseudocolour plot of product concentration relative to maximum reactant concen-
tration at t = 200s. Flow direction is from left to right.
D R A F T
May 31, 2012, 12:44pm
D R A F T
1
2
3
Simulation of bimolecular and trimolecular reactive transport in porous media using colocation probability functions. 1
1
J. M. Barnard, C. E. Augarde,
C. E. Augarde, School of Engineering and Computing Sciences, Durham University, South Road, Durham, DH1 3LE, UK. ([email protected]) J. M. Barnard, School of Engineering and Computing Sciences, Durham University, South Road, Durham, DH1 3LE, UK. ([email protected]) 1
School of Engineering and Computing
Sciences, Durham University, South Road, Durham, DH1 3LE, UK.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 2 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 4
Abstract.
Random walk particle tracking methods have shown them-
5
selves to be effective for describing mass transport in porous media due to
6
their ability to completely remove the problem of numerical dispersion. Sim-
7
ulation of reactions between different chemical species is a much less advanced
8
topic than in other methods of simulating mass transport in porous media
9
and has so far been restricted to simple reactions of the form A + B → C.
10
Here a method for simulating reactions with a stochastic particle tracking
11
framework is developed and discussed with respect to previous methods. This
12
method is based on colocation probability functions and is capable of deter-
13
mining the likelihood of the colocation of any number of particles with dif-
14
ferent arbitrary, Gaussian location probability functions. The method is val-
15
idated against experimental data from Gramling et al. [2002] and then used
16
in an expanded A + B + C → P form as a proof of concept for larger, more
17
complicated systems.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X-3
1. Introduction 18
The accurate simulation of reactive mass transport through porous media is an im-
19
portant capability for the study of various different environmental phenomena such as
20
contaminant transport within soil. As part of the multi-disciplinary, ROBUST project
21
(Regeneration Of Brownfield Using Sustainable Technology), which aims to develop a new
22
remediation technologies for contaminated land a new computational model for reactive
23
mass transport in porous media has been developed using a particle tracking based method
24
for simulating transport in tandem with a finite element porous media flow model.
25
26
The physical basis for porous media transport models is generally the advectiondispersion-reaction equation (ADRE) [Bear and Cheng, 2010];
27
∂C + v · ∇C = ∇ · (D∇C) + Q ∂t
28
where C is solute concentration, t is time, v is pore water velocity, D is a combined
29
dispersion-diffusion parameter, Q represents all of the sources and sinks for chemical
30
species including reactions and ∇ is the differential operator. There are many studies
31
into reactive mass transport through porous media that use either analytical or numerical
32
solutions to the ADRE as the basis of their work [Berkowitz et al., 2006; Edery et al.,
33
2010; Gramling et al., 2002].
(1)
34
Most models used to simulate porous media flow use a continuum approach to describe
35
the medium and its pore fluids. This works well for fluid flow problems but once solute
36
transport is added to the problem it becomes unviable for a range of reasons. Describing
37
the transport of a solute concentration field using a continuum based approach such as the
38
finite element method leads to smearing or numerical dispersion of the concentration field
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 4 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 39
parallel to the mean flow vector particularly in areas with high concentration gradients.
40
Though the extent of this effect can be reduced by refining the spatial resolution of the
41
simulation this is an inefficient and inelegant way to deal with the problem.
42
Another disadvantage of the continuum approach to describing mass transport in porous
43
media is that it cannot describe small scale variation in the solute concentration field. It
44
therefore assumes that reactants are mixed perfectly on the element scale. Experiment
45
[Gramling et al., 2002] and pore scale simulation [Cao and Kitanidis, 1998] have both
46
shown that the assumption of perfect mixing overpredicts the amount of a reactant pro-
47
duced.
48
Particle tracking techniques such as the discrete time random walk (DTRW) and con-
49
tinuous time random walk (CTRW) techniques completely eliminate the problem of nu-
50
merical dispersion from simulations of the ADRE and have been used in a wide range
51
of fields [Montroll and Weiss, 1965; Prickett et al., 1981]. Other advantages of particle
52
tracking methods are that they have been shown to be adept at describing anomalous
53
or non-Fickian transport [Berkowitz et al., 2006] and imperfect mixing of reactants in a
54
porous medium [Edery et al., 2010].
55
The application of the DTRW to the solution of the ADRE as set out by Prickett et al.
56
[1981] describes the concentration field of a chemical species using a field of particles.
57
These particles behave independently of all other particles of the same species and can be
58
located at any point within the domain.
59
In the DTRW method particles are advected according to the porous media flow field
60
and then walked a random distance ∆L over a fixed time step ∆t. The particles can take
61
any location within the domain unlike simulations using the CTRW method in which
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X-5
62
particle locations are limited to a lattice of locations. Simulating reactive mass transport
63
using either of the random walk methods is a more convoluted process than in any of the
64
continuum based approaches. As particles are treated as indivisable units of solute they
65
are reacted on a particle to particle basis. For example in the simulation of a reaction of
66
the form A + B → C, if a reaction takes place then an A and a B particle are destroyed
67
and a new C particle is formed.
2. Methods 68
In the random walk method each particle represents a ‘packet’ of the chemical species
69
rather than an individual particle, but is treated as though this is the case. The packet of a
70
chemical species discribed by a particle cannot change shape due to the effects of dispersion
71
or react partially with a packet of another species. These particles are advected in the
72
direction of the flow field and then dispersed a random distance parallel and perpendicular
73
to the direction of flow. Particle reactions are then performed on a particle by particle
74
basis.
75
The evolution of the particle field is performed in three stages; advection, dispersion
76
and reaction. The advection stage is performed by moving each particle according the
77
average flow field through the porous medium. This average flow field is calculated from
78
nodal pressure values using the methods of Cordes and Kinzelbach [1992] and Pollock
79
[1988]. These nodal pressure values could be determined using any continuum based
80
approach describing porous media flow. Here the finite element method is used with
81
bilinear quadrilateral elements.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 6 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
2.1. Discrete Time Random Walk Method 82
In the Pollock method particles are moved according to the flow field at their initial
83
location. The explanation of the method, from Pollock [1988], is for streamlines in two
84
dimensions. The expansion to three dimensions is trivial and so will not be given here.
85
The x and y components of a particle’s velocity are given as;
86
vxp =
(2)
vyp
(3)
87
1 (vx − vx1 )(xp − x1 ) + vx1 ∆x 2 1 = (vy − vy1 )(yp − y1 ) + vy1 ∆y 2
88
where vxp is the x component of the velocity of the flow field at the particle location,
89
∆x is the width of the element in which the particle is currently located, vx1 and vx2
90
are the x components of the flow field velocity at x1 and x2 , respectively, of the element,
91
xp is the x coordinate of the particle location and x1 and x2 are x coordinates of the
92
element boundaries perpendicular to x. The y components and coordinates of the above
93
parameters are labelled accordingly. This method was initially developed for use with
94
the finite difference method. It can be used unaltered with the finite element method if
95
regular bilinear quadrilateral elements are being used. If isoparametric elements are being
96
used then the global coordinates must be transformed so the element is square in the local
97
coordinate system.
98
The differential with respect to time of the x component of particle velocity is; dvx dt
99
!
= p
dvx dx
!
dx dt
!
.
(4)
p
100
The differential of vx with respect to x = (vx2 −vx1 )/∆x and the change of particle position
101
with respect to time is the definition of vxp . Therefore the simplification of (4) gives; dvx dt
102
D R A F T
!
= p
1 (vx − vx1 )vxp ∆x 2
May 31, 2012, 12:44pm
(5)
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 103
which can be rearranged to give; 1 1 dvxp = (vx − vx1 )dt. vxp ∆x 2
104
105
X-7
(6)
Integrating (6) between times t1 and t2 (t2 > t1 ) gives, "
#
106
vxp (t2 ) 1 ln = (vx − vx1 )∆t vxp (t2 ) ∆x 2
107
where ∆t = t2 − t1 . Rearrangement of (7) gives travel time from any x position within an
108
element to any other as; ∆tx =
109
110
∆x vx2 − vx1
!
"
(7)
#
vx2 ln . vx1
(8)
Analogous statements can be made for a particle’s movement with respect to y.
111
The first step in calculating the movement of a particle is to determine which element
112
edge it will leave through. This is done by calculating the travel times between a particle’s
113
location and each element edge using (8). Assuming the particle is not on an element
114
boundary the smallest positive travel time will correspond to the edge through which the
115
particle will exit. If the particle is on an element boundary then the travel time to that
116
boundary is obviously zero but the particle should only exit through that element edge if
117
the flow direction is towards it. To address this, if a particle is on an element boundary,
118
the element which it is in is determined by the direction of the pore fluid flux at that
119
boundary, the particle being placed in the element it is about to advect into rather than
120
advect out of. In addition if a particle is on an element edge, travel times of zero are
121
ignored when determining which edge the particle will exit through.
122
Dispersive transport is then simulated by the random walk of the particles. Each particle
123
is moved a random distance in the directions parallel and perpendicular to the mean flow
124
vector along which the particle was advected. A pseudo random number generator with D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 8 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 125
an N[0,1] distribution produced using a Box-M¨ uller transform on a U[-1,1] pseudo random
126
number generator is used. The generator is then scaled to the dispersivity and the time
127
step. The location probability function (LPF) for each particle is therefore Gaussian in
128
form with a standard deviation, σ, dependent on dispersivity and time step length. The
129
time dependence of dispersive transport can take one of two forms, either the Fickian form
130
where σ = (Dh ∆t)−0.5 or a truncated power law form (TPL) [Edery et al., 2010; Berkowitz
131
et al., 2006] where σ = e−t/t2 (1 + t/t1 )−1−β where β is the power law exponent, t1 is the
132
characteristic transition time and t2 is the cutoff time to Fickian transport. This TPL
133
form of time dependency was used with the CTRW method by Dentz et al. [2004], hence t1
134
being labelled in accordance with CTRW terminology. The TPL form of time dependency
135
has been used due to its ability to describe anomalous or non-Fickian transport which
136
occurs in heterogenous porous media [Berkowitz et al., 2000; Dentz et al., 2004; Edery
137
et al., 2010]. Edery et al. [2010] have shown that use of the TPL form of time dependency
138
can produce closer fits to experimental data for random walk based simulations.
139
Reactions in the random walk methods can be handled in several ways. Edery et al.
140
[2010] react particles if their separation is less than a reaction radius R. This reaction
141
radius is set empirically to fit the experimental data. Benson and Meerschaert [2008]
142
determine the likelihood that particles will be colocated at the end of the time step as a
143
function of current separation and time step size, or that for the next time step relative
144
motion will equal current separation. If this probability is greater than the value of a
145
random number from a U[0,1] generator then the particles are reacted. The probability of
146
this occurance is given by the convolution of two Gaussian LPFs which is also a Gaussian
147
function. Benson and Meerschaert [2008] approximate this Gaussian colocation probabil-
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X-9
148
ity function (CLPF) with a tent function to increase computational efficiency. This works
149
well when the location probability functions only have a single value for σ. If the LPFs
150
are of the more complicated form, such as when longitudinal and transverse dispersivity
151
are not equal, then determining the properties of the CLPF and the tent function to be
152
used become more complicated. This reduces any computational advantage gained from
153
the use of a tent function. 2.2. Particle Based Reaction Methods
154
In this study reaction is handled in a stochastic fashion based on the probability of one
155
particle of each of the reactant species being colocated over a time step ∆t and is similar
156
in some respects to the method used by Benson and Meerschaert [2008]. A particle’s
157
probability location function is of the same form as the distribution of the pseudo random
158
number generator used in the dispersive transport step. The probability of a pair of
159
particles meeting at a certain point will be the convolution of each of their LPFs. As the
160
LPFs have Gaussian forms their convolution is trivial. The complete definite integral of
161
this convolution or colocation probability function will then give the probability that the
162
particles will meet at any location during the time step. This method will initially be
163
demonstrated for the convolution of two arbitraty Gaussian functions in two dimensions.
164
Expansion to the convolution of any number of arbitrary Gaussian functions will then be
165
shown to be trivial.
166
167
The location probability function of a particle is expressed as an arbitrary two dimensional Gaussian function with its long axis parallel to the mean flow vector; 2 +2b(x−x
z = He−(a(x−x0 )
168
D R A F T
0 )(y−y0 )+c(y−y0 )
May 31, 2012, 12:44pm
2)
(9)
D R A F T
X - 10 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 169
where; a=
170
cos2 θ sin2 θ + 2σx2 2σy2
(10)
sin2θ cos2θ + 4σx2 4σy2
(11)
sin2 θ cos2 θ + 2σx2 2σy2
(12)
171
b=−
172
173
c=
174
175
and H is the maximum value of z located at (x0 , y0 ) and θ is the angle of the mean flow
176
vector from the y axis..
177
The convolution of two of these functions is given by;
178
z = H1 H2 exp(−(a1 (x − x1 )2 + a2 (x − x2 )2 2b1 (x − x1 )(y − y1 ) +2b2 (x − x2 )(y − y2 ) + c1 (y − y1 )2 + c2 (y − y2 )2 ))
179
180
which can be rearranged to give; 2 +by 2 +cxy+dx+f y+g
z = Heax
181
182
(13)
(14)
where;
183
a = (a1 + a2 ),
(15)
b = (c1 + c2 ),
(16)
c = 2(b1 + b2 ),
(17)
d = −2(a1 x1 + a2 x2 + b1 y1 + b2 y2 ),
(18)
f = −2(c1 y1 + c2 y2 + b1 x1 + b2 x2 ),
(19)
g = a1 x21 + a2 x22 + c1 y12 + c2 y22 + 2b1 x1 y1 + 2b2 x2 y2 ,
(20)
184
185 186
187 188
189 190
191 192
193
194
and
195
D R A F T
H = H1 H2 . May 31, 2012, 12:44pm
(21) D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 11
196
197
The definite integral of the above function is equal to; 1
2Hπ(4ab − c2 )− 2 e
198
−4abg+af 2 +bd2 +c2 g−cdf 4ab−c2
(22)
199
where, for the convolution of n arbitrary Gaussian functions each centred at (xn , yn ), the
200
coefficients a, b, c, d, f and g are given by; a=
201
202
b=
203
204
c=
205
n X
cos2 θi sin2 θi + 2 2σy2i i=1 2σxi
n X
sin2 θi cos2 θi + 2 2σy2i i=1 2σxi
n X i=1
206
d = −2
207
n X
"
i=1 208
f = −2
209
n X
"
i=1 210
211
g=
n X i=1
"
(23)
−
(24)
sin2θi sin2θi + 2σx2i 2σy2i
(25)
cos2 θi sin2 θi sin2θi sin2θi + xi + − + yi 2 2 2σxi 2σyi 4σx2i 4σy2i
(26)
sin2θi sin2θi sin2 θi cos2 θi + y + − + xi i 2σx2i 2σy2i 4σx2i 4σy2i
(27)
#
#
"
#
"
#
cos2 θi sin2 θi 2 sin2 θi cos2 θi 2 sin2θi sin2θi + x + + yi + 2 − + xi yi (28) i 2 2 2 2 2σxi 2σyi 2σxi 2σyi 4σx2i 4σy2i #
"
#
"
#
212
The maximum value of the CLPF, H, is set to one so that the maximum probability of
213
colocation is one and this occurs when the particles are colocated at the beginning of the
214
time step.
3. Results & Discussion 215
The model outlined above will be verified against experimental data from Gramling
216
et al. [2002] showing the transport and reaction of two solutes to form a thrid. The
217
experiment was performed using a box 36 cm long, 5.5 cm high and 1.8 cm deep filled
218
with cryolite sand. The non-reactive transport experiment starts with the sand saturated
219
with a 0.02 M solution of sodium ethylenediaminetetraacetic acid (NaEDTA2− ) which
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 12 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 220
is then displaced by a 0.02 M solution of CuSO4 at a flow rate of 2.67 mL/min, which
221
gives an average pore fluid velocity of 1.21 x 10−2 cm/s. The formation of the product
222
CuEDTA2− was then measured [Figure 1].
223
The experiment is simulated as a two dimensional domain 36 cm long and 5.5 cm high.
224
The flow field exists at a steady state producing a pore fluid velocity of 1.21 x 10−2 cm/s.
225
The left and right hand side boundaries are set as fixed pressure boundaries and top and
226
bottom boundaries are set as no flow boundaries. This initial concentration C 0 of each
227
reactant corresponds to a particle density of 225 per cm2 .
228
The simulation outputs correspond well to the experimental results from Gramling
229
et al. [2002] [Figures 2, 3 and 4]. The shape and peak concentration both match the
230
experimental data well. To produce these results the reaction probability has not been
231
empirically fitted to the experiment like the method used by Edery et al. [2010]. For a
232
Fortran 90 implementation of the model the run time to reach the situation shown in
233
Figure 4 with a time step length of 1 s on notebook computer with a 2.10 GHz, dual core
234
processor and 3 GB of RAM was 1505 seconds. Note that due to the probabilistic nature
235
of the model runtimes are variable. Run times are with 10% of those for simulations using
236
reaction radius method in an otherwise identical code.
237
The overall process of the reaction algorithm is; select particle A, determine which
238
B particles are within a set radius R of the selected A, loop over the set of B and for
239
each calculate (22) and compare to the output a U[0,1] pseudorandom number generator.
240
The checking of particles against one another to determine the probability of a reaction
241
occuring is the most computationally expensive part of the simulation and so certain
242
tactics were implemented within the code to attempt to minimise this. The set of B
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 13
243
less than R away from the selected A are placed in proximity order and so also order of
244
decreasing probability of reaction. In addition if the probability of a reaction occurring
245
ever drops below 10−3 then the loop over the set of B is terminated. Together these
246
additions to the code reduce the number of times that (22) is calculated. The checking of
247
each B particle to see if it is within R of the selected A is computationally expensive if the
248
entire set of B particles is checked. Instead the set of B particles is divided into subsets of
249
the particles within a set distance of each node at the beginning of the subroutine. The
250
closest node to the selected A is then determined and only the subset of B associated with
251
the same node are checked against.
252
Figures 5 & 6 show the variation of the solution at t = 916s with varying time step.
253
For both ∆t = 1s [Figure 3] and ∆t = 2s [Figure 6] the fit to the experiment data is
254
good, though at ∆t = 0.5s [Figure 5] the simulation slightly over-predicts the amount
255
of CuEDTA produced. The principal advantage of the CLPF method is that, for the
256
experiment simulated here, it does not need fitting empirically, unlike the reaction ra-
257
dius method. It scales automatically with the diffusivity and dispersivity of the situation
258
including situations with different logitudinal and transverse dispersivities, thereby au-
259
tomatically adjusting itself to variable flow fields. The stated advantage of both the
260
reaction radius and tent function methods, that they reduce run times by removing the
261
computationally expensive parts of the CLPF method has been shown to be minimal.
262
The second simulation run is used as a proof of concept for reaction involving a higher
263
number of reactants, in this case the reaction is of the form A + B + C → D. The basic
264
setup is the same as Figure 1 with a domain filled with a solution of one of the reactants
265
and the other two reactants flowing in at the same concentration. The setup is altered so
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 14 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS 266
that the domain is 20 cm high and 30 cm long and the inflow of reactants into the domain
267
is different. The other reactants then flow into the domain under a constant flow regime
268
from the left hand side. One reactant enters from the top fifth of the left hand side and
269
the third from the lowest fifth. The average pore fluid velocity is 1.21 x 10−1 cm/s and the
270
diffusivities of the reactants and products are all 1.75 x 10−2 . The initial concentration
271
of A in the domain is 900 particles per cm2 Figures 7 & 8 show concentrations of D at
272
100s and 200s respectively. The concentration field is a Gaussian surface weighted by
273
particle locations. Product concentrations relative to reaction concentrations are lower
274
than those in previous example due to the much lower probability of the three particles
275
being colocated than two.
276
Benson and Meerschaert [2008] state that for the sake of computational efficiency re-
277
actions of the type A + B + C → P should be handled as two consecutive reactions, A
278
+ B → D and C + D → P. Whilst it is true that the simulation of an A + B + C → P
279
reaction is computationally demanding it is not necessarily true than a different method
280
should be used. By splitting A + B + C → P into A + B → D and C + D → P the
281
eventuality of D particles existing unreacted in the domain becomes a possibility. If this
282
is possible in the chemical reaction being simulated then it is not a problem, but if the
283
reaction can only take in A,B and C and produce P then it becomes a problem. Also
284
the probability of three particles being colocated is much less than the probability of two
285
particles being colocated, reducing the reaction rate. Which reaction rate is the true one
286
would then affect which would be the better method to use.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 15
4. Conclusion 287
The colocation probability function method for simulating reaction in particle tracking
288
based models of reactive mass transport in porous media flow has been expanded from
289
the version presented by Benson and Meerschaert [2008]. This expansion has included
290
the capability to simulate reactions between any number of particles, each with different,
291
arbitrary Gaussian probability location functions. This method has also been tested
292
against experimental data from Gramling et al. [2002] where it appears to perform well.
293
The simulation does not require the variation of any empirical factors to ensure a good
294
fit to the experiment. As a proof of concept, the multicomponent reaction method has
295
been tested on a reaction of the form A + B + C → D. The overall purpose of this
296
work has been to develop a method for simulating reaction in particle tracking based
297
models of transport in porous media that is capable of simulating any of the reactions
298
that take place within the subsurface. The expansion of this method to include other
299
capabilities such as the simulation of reversible reactions, radioactive solutes, sorbtion of
300
solutes to the soil skeleton, catalysed and non-instantaneous reactions is trivial from a
301
programming perspective. Therefore a framework which is capable of describing a wide
302
range of chemical processes and so a large proportion of the processes involved within soil
303
chemistry has been developed.
304
305
Acknowledgments. This research is part of the Regenerating Brownfield Using Sustainable Technology project. EPSRC Challenging Engineering Grant: EP/G028958/1.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 16 BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
References 306
J. Bear and A. H.-D. Cheng. Modeling Groundwater Flow and Contaminant Transport,
307
volume 23 of Theory and Application of Transport in Porous Media. Springer, 2010.
308
D.A. Benson and M.M. Meerschaert. Simulation of chemical reaction via particle track-
309
ing: Diffusion-limited versus thermodynamic rate-limited regimes. Water Resources
310
Research, 44:W12201, 2008.
311
B. Berkowitz, H. Scher, and S.E. Silliman. Anomalous transport in laboratory-scale,
312
heterogeneous porous media. Water Resources Research, 36(1):149–158, 2000. ISSN
313
0043-1397.
314
B. Berkowitz, A. Cortis, M. Dentz, and H. Scher. Modeling non-Fickian transport in
315
geological formations as a continuous time random walk. Reviews of Geophysics, 44(2),
316
2006. ISSN 8755-1209.
317
318
J. Cao and P. K. Kitanidis. Pore-scale dilution of conservative solutes: An example. Water Resources Research, 34:1941–1949, 1998.
319
C. Cordes and W. Kinzelbach. Continuous groundwater velocity fields and path lines in
320
linear, bilinear, and trilinear finite elements. Water resources research, 28(11):2903–
321
2911, 1992. ISSN 0043-1397.
322
M. Dentz, A. Cortis, H. Scher, and B. Berkowitz. Time behavior of solute transport in
323
heterogeneous media: Transition from anomalous to normal transport. Advances in
324
Water Resources, 27(2):155–173, 2004.
325
326
Y. Edery, H. Scher, and B. Berkowitz. Particle tracking model of bimolecular reactive transport in porous media. Water Resources Research, 46(7):W07524, 2010.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 17
327
C.M. Gramling, C.F. Harvey, and L.C. Meigs. Reactive transport in porous media: A
328
comparison of model prediction with laboratory visualization. Environmental science
329
& technology, 36(11):2508–2514, 2002.
330
331
332
333
E.W. Montroll and G.H. Weiss. Random walks on lattices. II. Journal of Mathematical Physics, 6:167, 1965. D.W. Pollock. Semianalytical computation of path lines for finite-difference models. Ground Water, 26(6):743–750, 1988.
334
T.A. Prickett, T.G. Naymik, and C.G. Lonnquist. A’ Random Walk’ Solute Transport
335
Model for Selected Groundwater Quality Evaluations. Illinois State Water Survey Bul-
336
letin, 65, 1981.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 18
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
36 cm
0.02 M
NaEDTA2-
2.67 mL/min
Figure 1. i)
5.5 cm
0.02 M
CuSO4
Schematic diagram showing the experimental setup used by Gramling et al. [2002] t=n
ii)
t = n+1
and then simulated using the DTRW method. L
0.45
C
C
0.4
CuEDTA/C0
0.35
iii)
Simulation Experiment
x
0.3
t=n
iv)
t = n+1
0.25 0.2 0.15 0.1
C
x
L
C
0.05 0 0
x5
10
15
20
x
x (cm)
Figure 2. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 619 s.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
X - 19
0.45 Simulation Experiment
0.4
CuEDTA/C0
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0
5
10
15
20
x (cm)
Figure 3. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 916 s.
0.45 Simulation Experiment
0.4
CuEDTA/C0
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 5
10
15
20
25
x (cm)
Figure 4. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 1510 s.
D R A F T
May 31, 2012, 12:44pm
D R A F T
X - 20
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
0.45 Simulation Experiment
0.40
CuEDTA/C0
0.35 0.30
0.25 0.20 0.15
0.10 0.05 0.00 0
5
10
15
20
x (cm)
Figure 5. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 916 s for ∆t = 0.5s.
0.45 Series1 Series2
0.4
CuEDTA/C0
0.35 0.3
0.25 0.2 0.15
0.1 0.05 0 0
5
10
15
20
x (cm)
Figure 6. Comparsion of DTRW simulation and experimental data from Gramling et al. [2002] at t = 916 s for ∆t = 2s.
D R A F T
May 31, 2012, 12:44pm
D R A F T
BARNARD & AUGARDE: REACTIVE TRANSPORT USING COLOCATION PROBABILITY FUNCTIONS
Figure 7.
X - 21
Pseudocolour plot of product concentration relative to maximum reactant concen-
tration at t = 100s. Flow direction is from left to right.
Figure 8.
Pseudocolour plot of product concentration relative to maximum reactant concen-
tration at t = 200s. Flow direction is from left to right.
D R A F T
May 31, 2012, 12:44pm
D R A F T
