- Email: [email protected]
A grid-free particle tracking simulation for tracer dispersion in porous reservoir model
Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
Contents lists available at ScienceDirect
Journal of Unconventional Oil and Gas Resources journal homepage: www.elsevier.com/locate/juogr
A grid-free particle tracking simulation for tracer dispersion in porous reservoir model Arif Widiatmojo ⇑, Kyuro Sasaki, Amin Yousefi-Sahzabi, Ronald Nguele, Yuichi Sugai, Atsushi Maeda Department of Earth Resources Engineering, Kyushu University, Fukuoka, Japan
a r t i c l e
i n f o
Article history: Received 23 January 2015 Revised 4 May 2015 Accepted 25 May 2015 Available online 9 June 2015 Keywords: Tracer test Porous media Dispersion Numerical simulation Particle tracking method
a b s t r a c t Tracer test is a useful method to investigate various phenomena in geological porous media including groundwater contaminant transport, sweep efficiency and retention time in oil reservoir, reservoir characterization, fractures orientation assessment, as well as geothermal reservoir evaluation. Numerical methods are powerful tools in interpreting tracer test results. However, they are limited by computational restrictions which include finer grid requirements and small calculation steps. In this study, an analog model of a quarter five-spot porous reservoir was simulated by using random walk particle tracking method. This scheme used ‘method of images’ with pairs of injector–producer potential flow to generate the velocity vectors instead of conventionally solving Darcy’s equation to obtain grid velocities. Simulated breakthrough concentration profiles and flow visualization were compared with both experimental results and Eulerian-grid based finite volume simulation. The predicted breakthrough curves of tracer concentration were found to agree with experimental data sets. In addition to be free from numerical errors as often encountered in grid-based simulation, the proposed particle tracking model showed a faster computational time. Unlike the conventional grid method, this technique provides inherently smooth and continuous flow field at arbitrary position within the reservoir model. Ó 2015 Elsevier Ltd. All rights reserved.
Introduction The inter-well tracer test has been used to investigate the reservoir characteristics related to the effects of reservoir heterogeneity, sweep efficiency and convection–diffusion behavior of fluids. A typical tracer test is carried out by pulse (slug) release which consists of injecting a specific tracer fluid and monitoring its breakthrough profile at producer well. Throughout the flow from injector to producer, tracer fluid encounters various patterns inherent to reservoir characteristics. Thus, the breakthrough curves are considered to represent the porous characteristics of the reservoir within which the tracer test was conducted. Numerical simulations are essential analytical tools that provide means to interpret and analyze tracer test results. The conventional applied Eulerian-based transport models suffer from numerical dispersion and artificial oscillations, which are emphasized when advection is dominating. To overcome this problem, grid Peclet number must be chosen to be small enough to achieve numerical stability, resulting in long computational time (Salamon et al., 2006). Other methods have been developed as
⇑ Corresponding author at: 744 Motooka, Nishi-ku, Fukuoka-shi 819-0395, Japan. E-mail address: [email protected] (A. Widiatmojo). http://dx.doi.org/10.1016/j.juogr.2015.05.005 2213-3976/Ó 2015 Elsevier Ltd. All rights reserved.
alternative approach to aforementioned techniques, including Total Variation Diminishing (TVD) (Cox and Nishikawa, 1991), Method of Characteristic (Chiang et al., 1989) and Dynamic Mesh Refinement (Wolfsberg and Freyberg, 1994). On the other hand, Lagrangian-based particle tracking method treats species of solute mass as a large number of discrete particles. Not only requires no direct calculation of mass concentration from the discretized advection–dispersion equation, this method also inherently satisfies the mass conservation. The random walk particle tracking has been applied for the analysis of dispersion phenomena in porous media (Scheidegger, 1954; Josselin de Jong, 1958), based on the principal statistic of random Brownian movement (Einstein, 1905). Kinzelbach (1987, 1988) has noted unrealistic accumulation of particles in low permeability (i.e., low velocity) which may lead into significant concentration estimation error. Thus, correction terms were to be considered especially in the flow with high velocity variations or flow with strong injector and producer wells. Crane and Blunt (1999) and Obi and Blunt (2004) proposed the streamline based calculation of advection–dispersion in porous media which consisted of solving advection dispersion along the streamlines by using operator splitting technique. Chevalier and Banton (1999a,b) proposed the application of random walk method to the thermal energy storage in porous and fractured aquifers. In
76
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
Fig. 1. Top: illustration depicting reservoir model configuration; bottom: photograph of reservoir model.
their works, stored thermal energy are represented by particles which are moved within fractured aquifers. Practical applications of the dispersion in porous media may also involve kinetic reactions, radioactive decay or reactive chemical transport from two or more reactants which form new chemical products (Kinzelbach, 1988; Tompson, 1993; Cui et al., 2014). Cases by which the reactants are moving fluids or occurrence of chemical reactions between fluids and their surrounding matrix were reported. The applications of particle tracking method for the first case were primarily presented, for example by Cui et al. (2014). The reactants were treated as ensemble of moving particles in a flow field and their locations were recorded in order to obtain grid concentrations which were later used to calculate the reactive terms. Benson and Meerschaert (2008), Ding et al. (2013), and Paster et al. (2013, 2014) demonstrated the applicability of particle tracking method for bimolecular reaction by giving a probability value to different particles to react. It was able to perform direct calculation of reactions with no requirement of re-mapping particles’ into underlying grids to obtain local concentration. In this study, we examined the dispersion of tracer solute in a homogeneous-constrained quarter of a five spot system by constructing analog laboratory models. Two outcomes were taken into account from this analog model. First, the breakthrough curves of
arriving tracer concentration which undergo hydrodynamic dispersion through non-uniform velocity fields, and second, the spatial–temporal visualization of dispersed tracer. We compared the experimental results with numerical simulation by using particle tracking and advection–diffusion equation which was solved in a conventional finite volume. A quarter five spot model is characterized by its velocities’ non-uniformity, where velocities near the injector and producer wells appear to be high but are very low at two adjacent corners (stagnation zone). One of the most important steps of modeling solute transport in non-uniform flow fields is the generation of velocity field. For the conventional grid method, velocity vectors, u (m/s), are calculated within cell faces by solving Darcy’s equation:
u¼
k rp /l
ð1Þ
where k (m2) is permeability, / (–) is effective porosity of the porous system, l (Pa s) is the viscosity and p (Pa) are pressure. Furthermore, the particle tracking requires smooth-continuous velocity vectors. When this method is applied for particle tracking, velocity between cells must be interpolated (Crane and Blunt, 1999;
77
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
Obi and Blunt, 2004; Pollock, 1988). Moreover, errors may occur if the grid size is not small enough while velocity variations is high. In this paper, we proposed the superposition of sink–source potential flows (method of images) which provides inherently smooth and continuous velocity field. Moreover, not only this method eliminates possible errors due to velocities’ discontinuity at grid interfaces, but also provides faster computational time. Analog laboratory model The analog laboratory scale reservoir model with the dimension of L = 0.47 m, W = 0.30 m and T = 0.01 m was constructed from two acrylic plates (0.55 0.38 m) and sealed with rubber (0.01 m thickness) in between, resulting in reservoir dimension of 0.47 0.30 0.01 m. This reservoir which was filled with compacted glass beads to ensure perfect packing, was laid horizontally. After compaction, the model was weighed to obtain empty weight and then filled with a known volume of water until saturation was reached. The total weight with water was obtained by re-weighing the model after saturation. Based on this measurement, the porosity of the model was able to be estimated. Two holes were used as injector and producer for the solute as shown in Fig. 1. Sodium Chloride (NaCl)–water solute was used as tracer fluid and its concentration at the producer was monitored by measuring its conductivity. Prior injecting the tracer, a continuous water flooding was performed until saturation and steady flow was attained. Sodium Chloride was injected in very short period at injector to provide slug (instantaneous) tracer release. The conductivity of arriving solute at producer well was measured by using a conductivity meter. Visualization of dispersed tracer was performed by mixing a red-colored non-reactive agent with NaCl. The experiments considered two different flow rates, 20 ml/min and 25 ml/min.
where Dm is molecular diffusion coefficient (m2/s), I is identity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi matrix, juj ¼ u2x þ u2y is magnitude of velocity vector, aL and aT are longitudinal and transverse dispersivity respectively, ux and uy are velocity component in x and y direction respectively. In this paper, the velocity vectors are generated by two methods; the superpositioning of sink–source potential flow for particle tracking method and solving steady-state Darcy equation for the conventional finite volume scheme. The numerical scheme of conventional finite volume method was performed as (i) to generate the grids and assign flow properties, (ii) to solve the Darcy equation for each grid, (iii) to discritize convection–diffusion equation into algebraic equation and selection of discretization scheme and (iv) to solve the convection diffusion equation. In this study, the finite volume method as suggested by Patankar (1980) was used. The calculation domain of 0.47 m 0.30 m was discretized into square grid cell. The velocity vector at each grid cell interface is computed by using Darcy’s law and by giving the head gradient between injector and producer well and the permeability. From the various numerical differencing scheme, the second order upwind scheme was used for the present study. After several calculations, the 940 600 number of grids was selected as the optimum value to avoid numerical dispersion error. Initially (t = 0), the concentration at all control volumes other than injector are zero, while initial concentration, C0 is given to the injector well. The initial and boundary conditions are given as follow:
Cð0; 0; 0Þ ¼ C 0 ; Cðx; y; 0Þ ¼ 0 @C @x
ð0; y; tÞ ¼ @C ðL; y; tÞ ¼ @C ðx; 0; tÞ ¼ @C ðx; W; tÞ ¼ 0 @x @y @y
) ð4Þ
The no-slip boundary condition is applied for the velocity at the wall boundaries.
Numerical model
The particle tracking method
The transport of solute in porous media for steadyincompressible flow without adsorption or reaction of tracer in the reservoir formation can be described by the following advection–dispersion equation:
The particle tracking method takes advantage of the similarity between advection–diffusion equation (Eq. (1)) and the Fokker– Plank equation. This method simulates mass transport by discretizing mass concentration into numbers of particles. The movement of each particle is driven by advection due to drift vector and diffusion by random Brownian movement. The particle tracking scheme is represented by the following equation (Kinzelbach, 1987, 1988, 1990; Tompson and Gelhar, 1990; Kitanindis, 1994; Lichtner et al., 2002; Labolle et al., 1996):
@C þ r ðuCÞ ¼ r ðDrCÞ @t
ð2Þ
where C is concentration, t (s) is time, D is velocity dependent dispersion tensor, generally defined as:
D ¼ ðaT juj þ Dm ÞI þ ðaL aT Þ
uuT juj
ð3Þ
y1
(x1,0)
(x2,0)
Superimposed source (method of images)
Solid boundary (no-flow condition)
Fig. 2. Method of images applied to identical sources and method of generating solid (impervious) boundary.
pffiffiffiffiffiffi Xðt þ DtÞ ¼ XðtÞ þ ðu þ r DÞDt þ B Z Dt
ð5Þ
where X is particle location in Cartesian coordinate. The second term in right-hand side is the drift vector component of velocity vector and correction of dispersion gradient. The last term is the stochastic dispersion term with B is a displacement tensor equivalent to 2D = BBT and Z is vector of independent, normally distributed random number with zero mean and unit variance. The solution converges to true solution when number of particles is large enough to minimize oscillations. Since this particle tracking scheme requires no grid generation, the numerical restrictions encountered in conventional numerical method, i.e. grid Peclet number could be avoided. Initially, N0 number of non-reacting particles are released within a certain radius from injector location. Each particle is then moved by the convection due to fluid velocity and the dispersion. The convective movements for the next step are calculated from the velocity where the particles are located by using methods later described. The calculation is stopped when all particles arrived at producer well or when specified maximum number of calculation steps is reached. A particle is considered to reach at producer well
78
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
y axis
Producer (Sink)
Injector (Source)
x axis
Fig. 3. The computational domain for particle tracking method and the pattern of image well.
Fig. 4. Snapshots of eight different particle trajectories from injection (bottom-left corner) to the producer (upper-right corner); the magenta arrows are the velocity vectors generated from the method of images. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Table 1 Normalized CPU time required for calculation (finite volume method as a base comparison).
Normalized CPU time
Grid method (finite volume)
Particle tracking (N0 = 50,000; Np = 153)
Particle tracking (N0 = 15,000; Np = 91)
1 (reference)
0.49
0.29
and the calculation step for corresponding particle is terminated if its position is less than a specified radius from producer location. In order to match with the experimental breakthrough curve for each flow rate, the value of dispersion coefficients were adjusted by changing the dispersivities. The value of porosity was / = 0.38 as measured and similar to those commonly applied for close spheres packing (Bear, 1988). The boundary at wall was considered to be zero-flux boundary at which particles are bouncing back with the angle of reflection equal to the angle of incidence (Szymczak and Ladd, 2003). The reflection was assumed as fully elastic rebound at which the total distance is preserved.
Streamlines’ continuous velocity fields: the method of images The mathematical expression for continuity-conserving stream functions, w, at a point (x, y), in the flow field with a pair of source and sink located at (x0, y0) and (xi, yi) respectively is expressed as:
w¼
b y y0 b y yi tan1 tan1 2p 2p x x0 x xi
ð6Þ
The gradient at any point of a streamline defines the velocity vector at that point and its direction is perpendicular to potential lines. The method of images is a method by which streamlines, representing solid boundary, can be defined in a systematic arrangement (Bear, 1988). This scheme lies upon the fact that there is no existence of a flow across the separation between two or more identical sources or sinks. Consider two identical sources are placed at x axis, (x1, 0) and (x2, 0), as shown in Fig. 2, where x1 = x2, due to symmetry, the x-direction velocity component anywhere along y-axis must be zero. This can be used to treat the solid wall boundary. To create solid (impermeable) boundary of a quarter five-spot model, number of sink–source (injector–producer) pairs were arranged in a
79
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81 0.003
3.84E-04 3.80E-04
15 pairs
0.0025
3.76E-04
N/(N0Δτ)
3.72E-04 45 pairs
3.68E-04 153 pairs
N0=15,000 N0=50,000 0.0015
0.001
325 pairs
276 pairs
3.64E-04
0.0005
3.60E-04 0
50
100
150
200
250
300
350
Number of Pairs
0 0
Fig. 5. Velocity at (0.235, 0.15) for different number of pairs, Np = 15, 45, 153, 276, and 325.
systematic arrangement as shown in Fig. 3. The velocity at position (x, y) within flow field can be calculated by following equation by considering NP pairs of sources and sinks:
ux ¼ @w ¼ 2bp @y
Np X m¼1
uy ¼ @w ¼ 2bp @x
ðxx0 m Þ ðxx0 m Þ2 þðyy0 m Þ2
Np X m¼1
ðxx
ðyy0 m Þ ðxx0 m Þ2 þðyy0 m Þ2
ðxxi m Þ 2 2 i m Þ þðyyi m Þ
ðyyi m Þ 2 2 i m Þ þðyyi m Þ
ðxx
9 > > > > > =
ð7Þ
> > > > > ;
where w is stream function, Np is number of pairs, b (m2/s) is the source or sink strength, (x0, y0) and (xi, yi) are source and sink position respectively. It should be noted that Eq. (7) is the superposition of sources and sinks from Eq. (6). The strength of source/sink is evaluated by the following equation:
4Q /d
b¼
N0=4,000
0.002
ð8Þ
where Q (m3/s) is the injection flow rate, d (m) is the thickness of the reservoir model (0.01 m). The constant factor, 4, in Eq. (8), arises because of the fact that only a quarter of the flow domain was considered. Pairs of source and sink are positioned in such pattern to form the mirror wells of the injector and producer. Define a is the group number of pairs (a = 1, 2. . .), the minimum number of pairs to form a symmetrical arrangement of mirror well is 6 pairs (for a = 1). The dashed parallelograms in Fig. 3 show the arrangements of pairs with the smallest parallelogram (inner most) is 6 pairs system (a = 1). Prior applying whole numerical scheme, we provided an algorithm that generates those well images automatically
500
1000
1500
2000
2500
3000
Elapsed Time (s) Fig. 7. Tracer breakthrough curves from different number of particles (Q = 25 ml/ min, Np = 153 pairs).
according to number of pairs (a number) and the distance between injector and producer well. Additionally, the method of images imposes a non-zero velocity at the wall boundary due to the well images position. Hence, in order to satisfy the requirements of zero fluid velocity at fluid– solid wall boundary, we defined the flow of both velocity components normal and parallel to the wall to be zero, (ux = uy = 0 at x = 0, x = 0.47, y = 0 and y = 0.30). In this study, the fluid viscous effect at boundary layer adjacent to the wall has been neglected. Model verification The trajectories, obtained from 8 random particles from injector to producer, are shown in Fig. 4 along with the velocity vectors generated by the method of images. The small fluctuations as seen from each trajectory were found to be the result of the dispersion steps. Table 1 shows the CPU time required to perform calculations presented in the normalized value against CPU time for the grid method. The particle tracking method, using N0 = 50,000 and Np = 153 requires about half of the time required by grid method. Reduction of number of particles and pairs to N0 = 15,000 and Np = 91 pairs has further reduced the calculation time. We performed several sensitivity analyses on the numerical parameter used in the particle tracking. These analyses considered different values of a parameter at each set of simulations while other parameters were kept constant. Fig. 5 depicts the velocities at (0.235; 0.15) for 15, 45, 153, 276 and 325 pairs. Both curves are approaching asymptotic values as 0.002
0.003 15 pairs 45 pairs 66 pairs 153 pairs 276 pairs
N/(N0Δτ)
0.002
Δt=0.25s Δt =1s Δt =8s
0.0015
N/(N0 Δτ)
0.0025
0.0015
0.001
0.001 0.0005
0.0005
0
0
0
500
1000
1500
2000
2500
3000
Elapsed Time (s) Fig. 6. Tracer breakthrough curves from different number of pairs (Q = 25 ml/min, N0 = 20,000, dt = 0.5).
0
500
1000
1500
2000
2500
3000
3500
4000
Elapsed Time (s) Fig. 8. Tracer breakthrough curves from different time steps (Q = 20 ml/min, N0 = 20,000, Np = 153 pairs).
80
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
considering 276, 153, 66, 45 and 15 pairs. The breakthrough curves for particle tracking simulation were presented in the normalized value, N/(N0Ds) (s1), where Ds (s) is the counting time interval. The value of Ds is defined as:
0.002
Experiment 0.0015
N/(N0Δτ)
C/M
Grid Method Particle Tracking
Ds ¼ 0.001
0.0005
0 0
500
1000
1500
2000
2500
3000
3500
4000
Elapsed Time (s) Fig. 9. Comparison between particle tracking method with grid method (finite volume) and experimental results (20 ml/min, Np = 153 pairs, N0 = 40,000).
0.0025
0.002
N/(N0Δτ)
C/M
Experiment Grid Method Particle Tracking
0.0015
0.001
0.0005
0 0
500
1000
1500
2000
2500
3000
Elapsed Time (s) Fig. 10. Comparison between particle tracking method with grid method (finite volume) and experimental results (25 ml/min, Np = 153 pairs, N0 = 40,000).
the number of pairs increase. It suggests that more number of pairs will not cause any major velocity differences and the value is likely converged. Fig. 6 illustrates the particle tracking method
1 Dt 6 1 Dt Dt > 1
ð9Þ
Concentration breakthroughs reveal earlier arriving time as the number of pairs are reduced. However, this effect is diminished when larger number of pairs is considered, as can be seen from the cases of Np = 153 pairs and Np = 276 pairs. It is shown that the reduction of pairs to 45 has less significant effect on the curves even the earlier arriving tracers are noticeable at the front. However, a significant deviation can be seen for the case of 15 pairs. From these observations, it is clear that the number of pairs determines the accuracy of velocity prediction. The larger number of pairs will minimize the velocity error due to the geometrical position of well images. For further simulations, we used Np = 153 as optimum value for which the generated error was acceptably small and could be neglected. The concentration breakthrough curves obtained by considering 25 ml/min flow rate, Np = 153 pairs and different number of particles, N0 = 4,000; N0 = 15,000 and N0 = 50,000, as shown in Fig. 7 highlight relationship between number of particles and the smoothness of the breakthrough curves. The results in Fig. 7 suggest that the larger number of particles eliminates the effects of random fluctuations on the computed concentrations. This is particularly important if the calculation of local solute concentration is necessary when concentration-dependent chemical process is considered; for example chemical sorption or reaction of more than two chemical species. Fig. 8 shows the comparison of concentration breakthrough curves with different time steps for 20 ml/min flow rate. A close look at the curves Dt = 0.25 s and Dt = 8 s, suggests that the curve with Dt = 8 s has a slightly earlier arriving time and lesser concentration trailing compared with Dt = 0.25 s, which is believed to be caused by the longer ‘jump’ of a particle for each step as large time step is chosen. Figs. 9 and 10 illustrate the comparison of simulation results by particle tracking and conventional finite volume method with experimental results for 20 ml/min and 25 ml/min respectively.
5 min
10 min
15 min
Experimental results
Numerical calculation (grid model)
Particle Tracking
Fig. 11. Comparison of solute clouds between experimental results, finite volume and particle tracking method (Q = 20 ml/min).
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
The experimental data were presented in the normalized value C/M (s1), where M is defined as:
M¼
Z
T
C dt
ð10Þ
0
where T is total measuring time of arriving concentration. The simulations were done by evaluating longitudinal dispersivity aL = 1.2 103 m and transverse dispersivity, aT = aL/100 = 1.2 105 m for both 20 ml/min and 25 ml/min measurements. For 20 ml/min, the simulation results indicate some disagreements with the experimental results, although a consistency in both simulations were observed. Additionally, for 20 ml/min, it was found that numerical simulation could not reach a longer trailing concentration as indicated by the experimental results. The visualization of solute movements, carried out at 5 min, 10 min and 15 min after injection is presented in Fig. 11. Visualizations for both grid model and particle tracking revealed a good agreement with experimental visualizations. Instead of showing concentration by gradual color as in grid based simulation, the visualization of particle tracking was shown by plotting particles’ position. Hence, there are possibilities of overlapping which may lead into comparably ‘thicker’ profiles than experimental and grid based simulation results. It should be noted also that the ‘hook’ on the experimental results were caused by the viscous effect with the wall, whereas, it was not modeled in the numerical simulations as previously mentioned. It also explains the longer trailing concentration in experimental results compared to breakthrough concentrations obtained from numerical simulation, which was found to be greater for smaller flow rate (20 ml/min). Conclusion In this study, we demonstrated the free-grid simulation of tracer dispersion in a quarter of five spot model by particle tracking method. The velocity vectors were generated by mirror images of injector and producer wells. Compared to the conventional grid-based simulation, this method was found to be free from errors caused by numerical dispersion occurring in the conventional method. However, some numerical parameters which may lead to subsequent errors were observed at the breakthrough curves of arriving particles. The smoothness effect, caused by random fluctuations of arriving particles, were minimized by selecting higher number of particles. The accuracy of the generated velocity vectors was found to be strongly altered by the number of considered pairs. The more the number of pairs, the less the errors imputed to geometrical arrangement of injector and producer well images. This method can reduce calculation time compared to the generic grid-based finite volume method. Acknowledgement Authors would like to acknowledge Japan Society for the Promotion of Science (JSPS) for the financial supports (P14379/26-04379) and KAKENHI(14F04379).
81
References Bear, J., 1988. Dynamics of Fluids in Porous Media. Dover, New York. Benson, D.A., Meerschaert, M.M., 2008. Simulation of chemical reaction via particle tracking: diffusion-limited versus thermodynamic rate-limited regimes. Water Resour. Res. 44 (12), W12201. Chevalier, S., Banton, O., 1999a. Modelling of heat transfer with the random walk method. Part 1. Application to thermal energy storage in porous aquifers. J. Hydrol. 222 (1–4), 129–139. Chevalier, S., Banton, O., 1999b. Modelling of heat transfer with the random walk method. Part 2. Application to thermal energy storage in fractured aquifers. J. Hydrol. 222 (1–4), 140–151. Chiang, C.Y., Wheeler, M.F., Bedient, P.B., 1989. Method for simulation of groundwater solute transport. Water Resour. Res. 25 (7), 1541–1549. Cox, R.A., Nishikawa, T., 1991. A new total variation diminishing scheme for the solution of advective-dominant solute transport. Water Resour. Res. 27 (10), 2645–2654. Crane, M.J., Blunt, M.J., 1999. Streamline-based simulation of solute transport. Water Resour. Res. 35 (10), 3061–3078. Cui, Z., Welty, C., Maxwell, R.M., 2014. Modeling nitrogen transport and transformation in aquifers using a particle-tracking approach. Comput. Geosci. 70, 1–14. Ding, D. et al., 2013. Modeling bimolecular reactions and transport in porous media via particle tracking. Adv. Water Resour. 53, 56–65. Einstein, A., 1905. On the movement of small particles suspended in stationary liquids required by molecular-kinetic theory of heat. Ann. d Phys. 17, 549–560. Josselin de Jong, G., 1958. Longitudinal and transverse diffusion in granular deposits. In: Schotting, R.J., van Duijn, H.C.J., Verruijt, A. (Eds.), Trans. Am. Geophys. Union, 39. pp. 67–74. Kinzelbach, W., 1987. Methods for the simulation of pollutant transport in groundwater, a model comparison. In: Proceeding: Solving Ground Water Problems with Models, Conference and Exposition, vol. 1. Denver, Colorado, USA, pp. 656–675. Kinzelbach, W., 1988. The random walk method in pollutant transport simulation. Groundwater Flow Qual. Model. 224, 227–245. Kinzelbach, W., 1990. Simulation of pollutant transport in groundwater with the random walk method. In: Proceedings of the Groundwater Monitoring and Management. Dresden, pp. 265–279. Kitanindis, P.K., 1994. Particle-tracking equations for the solution of the advectiondispersion equation with variable coefficients. Water Resour. Res. 30 (11), 3225–3227. Labolle, E.M., Fogg, G.E., Tompson, A.F.B., 1996. Random-Walk simulation of transport in heterogeneous porous media: local mass-conservation problem and implementation methods. Water Resour. Res. 32 (3), 583–593. Lichtner, P.C., Kelkar, S., Robinson, B., 2002. New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking. Water Resour. Res. 38 (8), 21-1–21-16. Obi, E.-O., Blunt, M.J., 2004. Streamline-based simulation of advective–dispersive solute transport. Adv. Water Resour. 27 (9), 913–924. Paster, A., Bolster, D., Benson, D.A., 2013. Particle tracking and the diffusion-reaction equation. Water Resour. Res. 49 (1), 1–6. Paster, A., Bolster, D., Benson, D.a., 2014. Connecting the dots: semi-analytical and random walk numerical solutions of the diffusion–reaction equation with stochastic initial conditions. J. Comput. Phys. 263, 91–112. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. CRC Press. Pollock, D.W., 1988. Semianalytical computation of path lines for finite-difference models. Groundwater 26 (6), 743–750. Salamon, P., Fernàndez-Garcia, D., Gómez-Hernández, J.J., 2006. A review and numerical assessment of the random walk particle tracking method. J. Contam. Hydrol. 87 (3–4), 277–305. Scheidegger, A.E., 1954. Statistical hydrodynamics in porous media. J. Appl. Phys. 25 (8), 994. Szymczak, P., Ladd, A.J.C., 2003. Boundary conditions for stochastic solutions of the convection–diffusion equation. Phys. Rev. E 68 (3), 036704. Tompson, A.F.B., 1993. Numerical simulation of chemical migration in physically and chemically heterogeneous porous media. Water Resour. Res. 29 (11), 3709– 3726. Tompson, A.F.B., Gelhar, L.W., 1990. Numerical simulation of solute transport in three-dimensional, randomly heterogeneous porous media. Water Resour. Res. 26 (10), 2541–2562. Wolfsberg, V., Freyberg, D.L., 1994. Efficient simulation of single species and multispecies transport in groundwater with local adaptive grid refinement. Water Resour. Res. 30 (11), 2979–2991.
Contents lists available at ScienceDirect
Journal of Unconventional Oil and Gas Resources journal homepage: www.elsevier.com/locate/juogr
A grid-free particle tracking simulation for tracer dispersion in porous reservoir model Arif Widiatmojo ⇑, Kyuro Sasaki, Amin Yousefi-Sahzabi, Ronald Nguele, Yuichi Sugai, Atsushi Maeda Department of Earth Resources Engineering, Kyushu University, Fukuoka, Japan
a r t i c l e
i n f o
Article history: Received 23 January 2015 Revised 4 May 2015 Accepted 25 May 2015 Available online 9 June 2015 Keywords: Tracer test Porous media Dispersion Numerical simulation Particle tracking method
a b s t r a c t Tracer test is a useful method to investigate various phenomena in geological porous media including groundwater contaminant transport, sweep efficiency and retention time in oil reservoir, reservoir characterization, fractures orientation assessment, as well as geothermal reservoir evaluation. Numerical methods are powerful tools in interpreting tracer test results. However, they are limited by computational restrictions which include finer grid requirements and small calculation steps. In this study, an analog model of a quarter five-spot porous reservoir was simulated by using random walk particle tracking method. This scheme used ‘method of images’ with pairs of injector–producer potential flow to generate the velocity vectors instead of conventionally solving Darcy’s equation to obtain grid velocities. Simulated breakthrough concentration profiles and flow visualization were compared with both experimental results and Eulerian-grid based finite volume simulation. The predicted breakthrough curves of tracer concentration were found to agree with experimental data sets. In addition to be free from numerical errors as often encountered in grid-based simulation, the proposed particle tracking model showed a faster computational time. Unlike the conventional grid method, this technique provides inherently smooth and continuous flow field at arbitrary position within the reservoir model. Ó 2015 Elsevier Ltd. All rights reserved.
Introduction The inter-well tracer test has been used to investigate the reservoir characteristics related to the effects of reservoir heterogeneity, sweep efficiency and convection–diffusion behavior of fluids. A typical tracer test is carried out by pulse (slug) release which consists of injecting a specific tracer fluid and monitoring its breakthrough profile at producer well. Throughout the flow from injector to producer, tracer fluid encounters various patterns inherent to reservoir characteristics. Thus, the breakthrough curves are considered to represent the porous characteristics of the reservoir within which the tracer test was conducted. Numerical simulations are essential analytical tools that provide means to interpret and analyze tracer test results. The conventional applied Eulerian-based transport models suffer from numerical dispersion and artificial oscillations, which are emphasized when advection is dominating. To overcome this problem, grid Peclet number must be chosen to be small enough to achieve numerical stability, resulting in long computational time (Salamon et al., 2006). Other methods have been developed as
⇑ Corresponding author at: 744 Motooka, Nishi-ku, Fukuoka-shi 819-0395, Japan. E-mail address: [email protected] (A. Widiatmojo). http://dx.doi.org/10.1016/j.juogr.2015.05.005 2213-3976/Ó 2015 Elsevier Ltd. All rights reserved.
alternative approach to aforementioned techniques, including Total Variation Diminishing (TVD) (Cox and Nishikawa, 1991), Method of Characteristic (Chiang et al., 1989) and Dynamic Mesh Refinement (Wolfsberg and Freyberg, 1994). On the other hand, Lagrangian-based particle tracking method treats species of solute mass as a large number of discrete particles. Not only requires no direct calculation of mass concentration from the discretized advection–dispersion equation, this method also inherently satisfies the mass conservation. The random walk particle tracking has been applied for the analysis of dispersion phenomena in porous media (Scheidegger, 1954; Josselin de Jong, 1958), based on the principal statistic of random Brownian movement (Einstein, 1905). Kinzelbach (1987, 1988) has noted unrealistic accumulation of particles in low permeability (i.e., low velocity) which may lead into significant concentration estimation error. Thus, correction terms were to be considered especially in the flow with high velocity variations or flow with strong injector and producer wells. Crane and Blunt (1999) and Obi and Blunt (2004) proposed the streamline based calculation of advection–dispersion in porous media which consisted of solving advection dispersion along the streamlines by using operator splitting technique. Chevalier and Banton (1999a,b) proposed the application of random walk method to the thermal energy storage in porous and fractured aquifers. In
76
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
Fig. 1. Top: illustration depicting reservoir model configuration; bottom: photograph of reservoir model.
their works, stored thermal energy are represented by particles which are moved within fractured aquifers. Practical applications of the dispersion in porous media may also involve kinetic reactions, radioactive decay or reactive chemical transport from two or more reactants which form new chemical products (Kinzelbach, 1988; Tompson, 1993; Cui et al., 2014). Cases by which the reactants are moving fluids or occurrence of chemical reactions between fluids and their surrounding matrix were reported. The applications of particle tracking method for the first case were primarily presented, for example by Cui et al. (2014). The reactants were treated as ensemble of moving particles in a flow field and their locations were recorded in order to obtain grid concentrations which were later used to calculate the reactive terms. Benson and Meerschaert (2008), Ding et al. (2013), and Paster et al. (2013, 2014) demonstrated the applicability of particle tracking method for bimolecular reaction by giving a probability value to different particles to react. It was able to perform direct calculation of reactions with no requirement of re-mapping particles’ into underlying grids to obtain local concentration. In this study, we examined the dispersion of tracer solute in a homogeneous-constrained quarter of a five spot system by constructing analog laboratory models. Two outcomes were taken into account from this analog model. First, the breakthrough curves of
arriving tracer concentration which undergo hydrodynamic dispersion through non-uniform velocity fields, and second, the spatial–temporal visualization of dispersed tracer. We compared the experimental results with numerical simulation by using particle tracking and advection–diffusion equation which was solved in a conventional finite volume. A quarter five spot model is characterized by its velocities’ non-uniformity, where velocities near the injector and producer wells appear to be high but are very low at two adjacent corners (stagnation zone). One of the most important steps of modeling solute transport in non-uniform flow fields is the generation of velocity field. For the conventional grid method, velocity vectors, u (m/s), are calculated within cell faces by solving Darcy’s equation:
u¼
k rp /l
ð1Þ
where k (m2) is permeability, / (–) is effective porosity of the porous system, l (Pa s) is the viscosity and p (Pa) are pressure. Furthermore, the particle tracking requires smooth-continuous velocity vectors. When this method is applied for particle tracking, velocity between cells must be interpolated (Crane and Blunt, 1999;
77
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
Obi and Blunt, 2004; Pollock, 1988). Moreover, errors may occur if the grid size is not small enough while velocity variations is high. In this paper, we proposed the superposition of sink–source potential flows (method of images) which provides inherently smooth and continuous velocity field. Moreover, not only this method eliminates possible errors due to velocities’ discontinuity at grid interfaces, but also provides faster computational time. Analog laboratory model The analog laboratory scale reservoir model with the dimension of L = 0.47 m, W = 0.30 m and T = 0.01 m was constructed from two acrylic plates (0.55 0.38 m) and sealed with rubber (0.01 m thickness) in between, resulting in reservoir dimension of 0.47 0.30 0.01 m. This reservoir which was filled with compacted glass beads to ensure perfect packing, was laid horizontally. After compaction, the model was weighed to obtain empty weight and then filled with a known volume of water until saturation was reached. The total weight with water was obtained by re-weighing the model after saturation. Based on this measurement, the porosity of the model was able to be estimated. Two holes were used as injector and producer for the solute as shown in Fig. 1. Sodium Chloride (NaCl)–water solute was used as tracer fluid and its concentration at the producer was monitored by measuring its conductivity. Prior injecting the tracer, a continuous water flooding was performed until saturation and steady flow was attained. Sodium Chloride was injected in very short period at injector to provide slug (instantaneous) tracer release. The conductivity of arriving solute at producer well was measured by using a conductivity meter. Visualization of dispersed tracer was performed by mixing a red-colored non-reactive agent with NaCl. The experiments considered two different flow rates, 20 ml/min and 25 ml/min.
where Dm is molecular diffusion coefficient (m2/s), I is identity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi matrix, juj ¼ u2x þ u2y is magnitude of velocity vector, aL and aT are longitudinal and transverse dispersivity respectively, ux and uy are velocity component in x and y direction respectively. In this paper, the velocity vectors are generated by two methods; the superpositioning of sink–source potential flow for particle tracking method and solving steady-state Darcy equation for the conventional finite volume scheme. The numerical scheme of conventional finite volume method was performed as (i) to generate the grids and assign flow properties, (ii) to solve the Darcy equation for each grid, (iii) to discritize convection–diffusion equation into algebraic equation and selection of discretization scheme and (iv) to solve the convection diffusion equation. In this study, the finite volume method as suggested by Patankar (1980) was used. The calculation domain of 0.47 m 0.30 m was discretized into square grid cell. The velocity vector at each grid cell interface is computed by using Darcy’s law and by giving the head gradient between injector and producer well and the permeability. From the various numerical differencing scheme, the second order upwind scheme was used for the present study. After several calculations, the 940 600 number of grids was selected as the optimum value to avoid numerical dispersion error. Initially (t = 0), the concentration at all control volumes other than injector are zero, while initial concentration, C0 is given to the injector well. The initial and boundary conditions are given as follow:
Cð0; 0; 0Þ ¼ C 0 ; Cðx; y; 0Þ ¼ 0 @C @x
ð0; y; tÞ ¼ @C ðL; y; tÞ ¼ @C ðx; 0; tÞ ¼ @C ðx; W; tÞ ¼ 0 @x @y @y
) ð4Þ
The no-slip boundary condition is applied for the velocity at the wall boundaries.
Numerical model
The particle tracking method
The transport of solute in porous media for steadyincompressible flow without adsorption or reaction of tracer in the reservoir formation can be described by the following advection–dispersion equation:
The particle tracking method takes advantage of the similarity between advection–diffusion equation (Eq. (1)) and the Fokker– Plank equation. This method simulates mass transport by discretizing mass concentration into numbers of particles. The movement of each particle is driven by advection due to drift vector and diffusion by random Brownian movement. The particle tracking scheme is represented by the following equation (Kinzelbach, 1987, 1988, 1990; Tompson and Gelhar, 1990; Kitanindis, 1994; Lichtner et al., 2002; Labolle et al., 1996):
@C þ r ðuCÞ ¼ r ðDrCÞ @t
ð2Þ
where C is concentration, t (s) is time, D is velocity dependent dispersion tensor, generally defined as:
D ¼ ðaT juj þ Dm ÞI þ ðaL aT Þ
uuT juj
ð3Þ
y1
(x1,0)
(x2,0)
Superimposed source (method of images)
Solid boundary (no-flow condition)
Fig. 2. Method of images applied to identical sources and method of generating solid (impervious) boundary.
pffiffiffiffiffiffi Xðt þ DtÞ ¼ XðtÞ þ ðu þ r DÞDt þ B Z Dt
ð5Þ
where X is particle location in Cartesian coordinate. The second term in right-hand side is the drift vector component of velocity vector and correction of dispersion gradient. The last term is the stochastic dispersion term with B is a displacement tensor equivalent to 2D = BBT and Z is vector of independent, normally distributed random number with zero mean and unit variance. The solution converges to true solution when number of particles is large enough to minimize oscillations. Since this particle tracking scheme requires no grid generation, the numerical restrictions encountered in conventional numerical method, i.e. grid Peclet number could be avoided. Initially, N0 number of non-reacting particles are released within a certain radius from injector location. Each particle is then moved by the convection due to fluid velocity and the dispersion. The convective movements for the next step are calculated from the velocity where the particles are located by using methods later described. The calculation is stopped when all particles arrived at producer well or when specified maximum number of calculation steps is reached. A particle is considered to reach at producer well
78
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
y axis
Producer (Sink)
Injector (Source)
x axis
Fig. 3. The computational domain for particle tracking method and the pattern of image well.
Fig. 4. Snapshots of eight different particle trajectories from injection (bottom-left corner) to the producer (upper-right corner); the magenta arrows are the velocity vectors generated from the method of images. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Table 1 Normalized CPU time required for calculation (finite volume method as a base comparison).
Normalized CPU time
Grid method (finite volume)
Particle tracking (N0 = 50,000; Np = 153)
Particle tracking (N0 = 15,000; Np = 91)
1 (reference)
0.49
0.29
and the calculation step for corresponding particle is terminated if its position is less than a specified radius from producer location. In order to match with the experimental breakthrough curve for each flow rate, the value of dispersion coefficients were adjusted by changing the dispersivities. The value of porosity was / = 0.38 as measured and similar to those commonly applied for close spheres packing (Bear, 1988). The boundary at wall was considered to be zero-flux boundary at which particles are bouncing back with the angle of reflection equal to the angle of incidence (Szymczak and Ladd, 2003). The reflection was assumed as fully elastic rebound at which the total distance is preserved.
Streamlines’ continuous velocity fields: the method of images The mathematical expression for continuity-conserving stream functions, w, at a point (x, y), in the flow field with a pair of source and sink located at (x0, y0) and (xi, yi) respectively is expressed as:
w¼
b y y0 b y yi tan1 tan1 2p 2p x x0 x xi
ð6Þ
The gradient at any point of a streamline defines the velocity vector at that point and its direction is perpendicular to potential lines. The method of images is a method by which streamlines, representing solid boundary, can be defined in a systematic arrangement (Bear, 1988). This scheme lies upon the fact that there is no existence of a flow across the separation between two or more identical sources or sinks. Consider two identical sources are placed at x axis, (x1, 0) and (x2, 0), as shown in Fig. 2, where x1 = x2, due to symmetry, the x-direction velocity component anywhere along y-axis must be zero. This can be used to treat the solid wall boundary. To create solid (impermeable) boundary of a quarter five-spot model, number of sink–source (injector–producer) pairs were arranged in a
79
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81 0.003
3.84E-04 3.80E-04
15 pairs
0.0025
3.76E-04
N/(N0Δτ)
3.72E-04 45 pairs
3.68E-04 153 pairs
N0=15,000 N0=50,000 0.0015
0.001
325 pairs
276 pairs
3.64E-04
0.0005
3.60E-04 0
50
100
150
200
250
300
350
Number of Pairs
0 0
Fig. 5. Velocity at (0.235, 0.15) for different number of pairs, Np = 15, 45, 153, 276, and 325.
systematic arrangement as shown in Fig. 3. The velocity at position (x, y) within flow field can be calculated by following equation by considering NP pairs of sources and sinks:
ux ¼ @w ¼ 2bp @y
Np X m¼1
uy ¼ @w ¼ 2bp @x
ðxx0 m Þ ðxx0 m Þ2 þðyy0 m Þ2
Np X m¼1
ðxx
ðyy0 m Þ ðxx0 m Þ2 þðyy0 m Þ2
ðxxi m Þ 2 2 i m Þ þðyyi m Þ
ðyyi m Þ 2 2 i m Þ þðyyi m Þ
ðxx
9 > > > > > =
ð7Þ
> > > > > ;
where w is stream function, Np is number of pairs, b (m2/s) is the source or sink strength, (x0, y0) and (xi, yi) are source and sink position respectively. It should be noted that Eq. (7) is the superposition of sources and sinks from Eq. (6). The strength of source/sink is evaluated by the following equation:
4Q /d
b¼
N0=4,000
0.002
ð8Þ
where Q (m3/s) is the injection flow rate, d (m) is the thickness of the reservoir model (0.01 m). The constant factor, 4, in Eq. (8), arises because of the fact that only a quarter of the flow domain was considered. Pairs of source and sink are positioned in such pattern to form the mirror wells of the injector and producer. Define a is the group number of pairs (a = 1, 2. . .), the minimum number of pairs to form a symmetrical arrangement of mirror well is 6 pairs (for a = 1). The dashed parallelograms in Fig. 3 show the arrangements of pairs with the smallest parallelogram (inner most) is 6 pairs system (a = 1). Prior applying whole numerical scheme, we provided an algorithm that generates those well images automatically
500
1000
1500
2000
2500
3000
Elapsed Time (s) Fig. 7. Tracer breakthrough curves from different number of particles (Q = 25 ml/ min, Np = 153 pairs).
according to number of pairs (a number) and the distance between injector and producer well. Additionally, the method of images imposes a non-zero velocity at the wall boundary due to the well images position. Hence, in order to satisfy the requirements of zero fluid velocity at fluid– solid wall boundary, we defined the flow of both velocity components normal and parallel to the wall to be zero, (ux = uy = 0 at x = 0, x = 0.47, y = 0 and y = 0.30). In this study, the fluid viscous effect at boundary layer adjacent to the wall has been neglected. Model verification The trajectories, obtained from 8 random particles from injector to producer, are shown in Fig. 4 along with the velocity vectors generated by the method of images. The small fluctuations as seen from each trajectory were found to be the result of the dispersion steps. Table 1 shows the CPU time required to perform calculations presented in the normalized value against CPU time for the grid method. The particle tracking method, using N0 = 50,000 and Np = 153 requires about half of the time required by grid method. Reduction of number of particles and pairs to N0 = 15,000 and Np = 91 pairs has further reduced the calculation time. We performed several sensitivity analyses on the numerical parameter used in the particle tracking. These analyses considered different values of a parameter at each set of simulations while other parameters were kept constant. Fig. 5 depicts the velocities at (0.235; 0.15) for 15, 45, 153, 276 and 325 pairs. Both curves are approaching asymptotic values as 0.002
0.003 15 pairs 45 pairs 66 pairs 153 pairs 276 pairs
N/(N0Δτ)
0.002
Δt=0.25s Δt =1s Δt =8s
0.0015
N/(N0 Δτ)
0.0025
0.0015
0.001
0.001 0.0005
0.0005
0
0
0
500
1000
1500
2000
2500
3000
Elapsed Time (s) Fig. 6. Tracer breakthrough curves from different number of pairs (Q = 25 ml/min, N0 = 20,000, dt = 0.5).
0
500
1000
1500
2000
2500
3000
3500
4000
Elapsed Time (s) Fig. 8. Tracer breakthrough curves from different time steps (Q = 20 ml/min, N0 = 20,000, Np = 153 pairs).
80
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
considering 276, 153, 66, 45 and 15 pairs. The breakthrough curves for particle tracking simulation were presented in the normalized value, N/(N0Ds) (s1), where Ds (s) is the counting time interval. The value of Ds is defined as:
0.002
Experiment 0.0015
N/(N0Δτ)
C/M
Grid Method Particle Tracking
Ds ¼ 0.001
0.0005
0 0
500
1000
1500
2000
2500
3000
3500
4000
Elapsed Time (s) Fig. 9. Comparison between particle tracking method with grid method (finite volume) and experimental results (20 ml/min, Np = 153 pairs, N0 = 40,000).
0.0025
0.002
N/(N0Δτ)
C/M
Experiment Grid Method Particle Tracking
0.0015
0.001
0.0005
0 0
500
1000
1500
2000
2500
3000
Elapsed Time (s) Fig. 10. Comparison between particle tracking method with grid method (finite volume) and experimental results (25 ml/min, Np = 153 pairs, N0 = 40,000).
the number of pairs increase. It suggests that more number of pairs will not cause any major velocity differences and the value is likely converged. Fig. 6 illustrates the particle tracking method
1 Dt 6 1 Dt Dt > 1
ð9Þ
Concentration breakthroughs reveal earlier arriving time as the number of pairs are reduced. However, this effect is diminished when larger number of pairs is considered, as can be seen from the cases of Np = 153 pairs and Np = 276 pairs. It is shown that the reduction of pairs to 45 has less significant effect on the curves even the earlier arriving tracers are noticeable at the front. However, a significant deviation can be seen for the case of 15 pairs. From these observations, it is clear that the number of pairs determines the accuracy of velocity prediction. The larger number of pairs will minimize the velocity error due to the geometrical position of well images. For further simulations, we used Np = 153 as optimum value for which the generated error was acceptably small and could be neglected. The concentration breakthrough curves obtained by considering 25 ml/min flow rate, Np = 153 pairs and different number of particles, N0 = 4,000; N0 = 15,000 and N0 = 50,000, as shown in Fig. 7 highlight relationship between number of particles and the smoothness of the breakthrough curves. The results in Fig. 7 suggest that the larger number of particles eliminates the effects of random fluctuations on the computed concentrations. This is particularly important if the calculation of local solute concentration is necessary when concentration-dependent chemical process is considered; for example chemical sorption or reaction of more than two chemical species. Fig. 8 shows the comparison of concentration breakthrough curves with different time steps for 20 ml/min flow rate. A close look at the curves Dt = 0.25 s and Dt = 8 s, suggests that the curve with Dt = 8 s has a slightly earlier arriving time and lesser concentration trailing compared with Dt = 0.25 s, which is believed to be caused by the longer ‘jump’ of a particle for each step as large time step is chosen. Figs. 9 and 10 illustrate the comparison of simulation results by particle tracking and conventional finite volume method with experimental results for 20 ml/min and 25 ml/min respectively.
5 min
10 min
15 min
Experimental results
Numerical calculation (grid model)
Particle Tracking
Fig. 11. Comparison of solute clouds between experimental results, finite volume and particle tracking method (Q = 20 ml/min).
A. Widiatmojo et al. / Journal of Unconventional Oil and Gas Resources 11 (2015) 75–81
The experimental data were presented in the normalized value C/M (s1), where M is defined as:
M¼
Z
T
C dt
ð10Þ
0
where T is total measuring time of arriving concentration. The simulations were done by evaluating longitudinal dispersivity aL = 1.2 103 m and transverse dispersivity, aT = aL/100 = 1.2 105 m for both 20 ml/min and 25 ml/min measurements. For 20 ml/min, the simulation results indicate some disagreements with the experimental results, although a consistency in both simulations were observed. Additionally, for 20 ml/min, it was found that numerical simulation could not reach a longer trailing concentration as indicated by the experimental results. The visualization of solute movements, carried out at 5 min, 10 min and 15 min after injection is presented in Fig. 11. Visualizations for both grid model and particle tracking revealed a good agreement with experimental visualizations. Instead of showing concentration by gradual color as in grid based simulation, the visualization of particle tracking was shown by plotting particles’ position. Hence, there are possibilities of overlapping which may lead into comparably ‘thicker’ profiles than experimental and grid based simulation results. It should be noted also that the ‘hook’ on the experimental results were caused by the viscous effect with the wall, whereas, it was not modeled in the numerical simulations as previously mentioned. It also explains the longer trailing concentration in experimental results compared to breakthrough concentrations obtained from numerical simulation, which was found to be greater for smaller flow rate (20 ml/min). Conclusion In this study, we demonstrated the free-grid simulation of tracer dispersion in a quarter of five spot model by particle tracking method. The velocity vectors were generated by mirror images of injector and producer wells. Compared to the conventional grid-based simulation, this method was found to be free from errors caused by numerical dispersion occurring in the conventional method. However, some numerical parameters which may lead to subsequent errors were observed at the breakthrough curves of arriving particles. The smoothness effect, caused by random fluctuations of arriving particles, were minimized by selecting higher number of particles. The accuracy of the generated velocity vectors was found to be strongly altered by the number of considered pairs. The more the number of pairs, the less the errors imputed to geometrical arrangement of injector and producer well images. This method can reduce calculation time compared to the generic grid-based finite volume method. Acknowledgement Authors would like to acknowledge Japan Society for the Promotion of Science (JSPS) for the financial supports (P14379/26-04379) and KAKENHI(14F04379).
81
References Bear, J., 1988. Dynamics of Fluids in Porous Media. Dover, New York. Benson, D.A., Meerschaert, M.M., 2008. Simulation of chemical reaction via particle tracking: diffusion-limited versus thermodynamic rate-limited regimes. Water Resour. Res. 44 (12), W12201. Chevalier, S., Banton, O., 1999a. Modelling of heat transfer with the random walk method. Part 1. Application to thermal energy storage in porous aquifers. J. Hydrol. 222 (1–4), 129–139. Chevalier, S., Banton, O., 1999b. Modelling of heat transfer with the random walk method. Part 2. Application to thermal energy storage in fractured aquifers. J. Hydrol. 222 (1–4), 140–151. Chiang, C.Y., Wheeler, M.F., Bedient, P.B., 1989. Method for simulation of groundwater solute transport. Water Resour. Res. 25 (7), 1541–1549. Cox, R.A., Nishikawa, T., 1991. A new total variation diminishing scheme for the solution of advective-dominant solute transport. Water Resour. Res. 27 (10), 2645–2654. Crane, M.J., Blunt, M.J., 1999. Streamline-based simulation of solute transport. Water Resour. Res. 35 (10), 3061–3078. Cui, Z., Welty, C., Maxwell, R.M., 2014. Modeling nitrogen transport and transformation in aquifers using a particle-tracking approach. Comput. Geosci. 70, 1–14. Ding, D. et al., 2013. Modeling bimolecular reactions and transport in porous media via particle tracking. Adv. Water Resour. 53, 56–65. Einstein, A., 1905. On the movement of small particles suspended in stationary liquids required by molecular-kinetic theory of heat. Ann. d Phys. 17, 549–560. Josselin de Jong, G., 1958. Longitudinal and transverse diffusion in granular deposits. In: Schotting, R.J., van Duijn, H.C.J., Verruijt, A. (Eds.), Trans. Am. Geophys. Union, 39. pp. 67–74. Kinzelbach, W., 1987. Methods for the simulation of pollutant transport in groundwater, a model comparison. In: Proceeding: Solving Ground Water Problems with Models, Conference and Exposition, vol. 1. Denver, Colorado, USA, pp. 656–675. Kinzelbach, W., 1988. The random walk method in pollutant transport simulation. Groundwater Flow Qual. Model. 224, 227–245. Kinzelbach, W., 1990. Simulation of pollutant transport in groundwater with the random walk method. In: Proceedings of the Groundwater Monitoring and Management. Dresden, pp. 265–279. Kitanindis, P.K., 1994. Particle-tracking equations for the solution of the advectiondispersion equation with variable coefficients. Water Resour. Res. 30 (11), 3225–3227. Labolle, E.M., Fogg, G.E., Tompson, A.F.B., 1996. Random-Walk simulation of transport in heterogeneous porous media: local mass-conservation problem and implementation methods. Water Resour. Res. 32 (3), 583–593. Lichtner, P.C., Kelkar, S., Robinson, B., 2002. New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking. Water Resour. Res. 38 (8), 21-1–21-16. Obi, E.-O., Blunt, M.J., 2004. Streamline-based simulation of advective–dispersive solute transport. Adv. Water Resour. 27 (9), 913–924. Paster, A., Bolster, D., Benson, D.A., 2013. Particle tracking and the diffusion-reaction equation. Water Resour. Res. 49 (1), 1–6. Paster, A., Bolster, D., Benson, D.a., 2014. Connecting the dots: semi-analytical and random walk numerical solutions of the diffusion–reaction equation with stochastic initial conditions. J. Comput. Phys. 263, 91–112. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. CRC Press. Pollock, D.W., 1988. Semianalytical computation of path lines for finite-difference models. Groundwater 26 (6), 743–750. Salamon, P., Fernàndez-Garcia, D., Gómez-Hernández, J.J., 2006. A review and numerical assessment of the random walk particle tracking method. J. Contam. Hydrol. 87 (3–4), 277–305. Scheidegger, A.E., 1954. Statistical hydrodynamics in porous media. J. Appl. Phys. 25 (8), 994. Szymczak, P., Ladd, A.J.C., 2003. Boundary conditions for stochastic solutions of the convection–diffusion equation. Phys. Rev. E 68 (3), 036704. Tompson, A.F.B., 1993. Numerical simulation of chemical migration in physically and chemically heterogeneous porous media. Water Resour. Res. 29 (11), 3709– 3726. Tompson, A.F.B., Gelhar, L.W., 1990. Numerical simulation of solute transport in three-dimensional, randomly heterogeneous porous media. Water Resour. Res. 26 (10), 2541–2562. Wolfsberg, V., Freyberg, D.L., 1994. Efficient simulation of single species and multispecies transport in groundwater with local adaptive grid refinement. Water Resour. Res. 30 (11), 2979–2991.