PARTICLE AND STREAMLINE NUMERICAL METHODS FOR CONSERVATIVE AND REACTIVE TRANSPORT SIMULATIONS IN POROUS MEDIA by PAULO ANDRES HERRERA RICCI M.Sc., University of Illinois at Urbana-Champaign, 2003 Ingeniero Civil, Universidad de Chile, 2001 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Geological Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) November 2009 cPaulo Andres Herrera Ricci, 2009 Abstract Reactive transport modeling has become an important tool to study and understand the transport and fate of solutes in the subsurface. However, the accurate simulation of reactive transport represents a formidable challenge because of the characteristics of flow, transport and chemical reactions that govern the migration of solutes in geological formations. In particular, solute transport in natural porous media is advection-controlled and disper- sion is higher in the direction of flow than in the transverse direction. Both characteristics create difficulties for traditional numerical schemes that result in numerical dispersion and/or spurious oscillations. While these errors can often be tolerated in conservative transport simulations, they can be amplified in presence of chemical reactions resulting in much larger errors or unstable solutions. In this thesis, new Lagrangian based methods to simulate conservative and reactive transport in porous media are investigated. First, the derivation of a new meshless approximation based on smoothed particle hydrodynamics (SPH) to simulate conserva- tive multidimensional solute transport, including advection and anisotropic dispersion, is presented. Second, a hybrid scheme that combines some of the advantages of streamline- based simulations and meshless methods and that allows simulating longitudinal and transverse dispersion without requiring a background grid is also derived. The numer- ical properties of both methods are analyzed analytical and numerically. Furthermore, both formulations are compared with existing numerical techniques in a set of two- and three-dimensional benchmark problems. It is demonstrated that the proposed schemes provide accurate and efficient solutions of physical transport processes in heterogeneous porous media and overcome most of the issues in existing numerical formulations. The new methods have the potential to remove or minimize numerical dispersion and grid orientation effects and, in the case ii of the hybrid streamline method, also eliminate spurious oscillations even in presence of large longitudinal to transverse dispersivity ratios. Therefore, the results presented in this thesis confirm that the Lagrangian formulations of solute transport investigated here are viable and compelling alternatives to simulate reactive transport versus more standard numerical techniques. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi Co-authorship Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Reactive Transport Modeling . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Conceptual Model . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.3 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.1 Particularities of Flow and Transport in Porous Media . . . . . 6 1.3.2 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.2.1 Mesh-based numerical methods . . . . . . . . . . . . . 10 1.3.2.2 Hybrid Eulerian-Lagrangian methods . . . . . . . . . . 12 1.3.2.3 Random walk particle tracking methods . . . . . . . . 12 1.3.3 Limitations of Current Numerical Methods . . . . . . . . . . . . 13 1.3.3.1 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.3.2 Monotonicity . . . . . . . . . . . . . . . . . . . . . . . 15 1.3.3.3 Performance . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4 Lagrangian Numerical Methods . . . . . . . . . . . . . . . . . . . . . . 17 1.4.1 Meshless Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4.2 Streamline-Based Simulations . . . . . . . . . . . . . . . . . . . 18 1.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 iv 1.6 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.7 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.8 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2 A Meshless Method to Simulate Solute Transport in Heterogeneous Porous Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.1.2 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 Monte Carlo SPH method . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.1 Time Integration . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.2.2 Accuracy and Spatial Resolution . . . . . . . . . . . . . . . . . 39 2.2.3 Mass Conservation . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.3 Numerical evaluation of the MC-SPH method . . . . . . . . . . . . . . 41 2.3.1 One-Dimensional Dispersion . . . . . . . . . . . . . . . . . . . . 41 2.3.2 Two-Dimensional Dispersion . . . . . . . . . . . . . . . . . . . . 42 2.3.2.1 Initial particle and concentration distribution . . . . . 46 2.3.2.2 Performance . . . . . . . . . . . . . . . . . . . . . . . . 48 2.3.2.3 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.3.3 Advection–Dispersion in Heterogeneous Porous Media . . . . . . 54 2.3.3.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.3.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.5 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3 Evaluation of Particle Approximations to Simulate Anisotropic Dis- persion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.2 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.3 Smoothed Particle Hydrodynamics (SPH) Approximation . . . . . . . . 75 3.3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.3.2 SPH Approximation for Tensorial Dispersion . . . . . . . . . . . 76 3.3.3 Monotonicity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.4 Particle Strength Exchange (PSE) Approximation . . . . . . . . . . . . 80 3.5 Numerical Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.5.1 Simulation Cases . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.5.2 Equispaced Particles . . . . . . . . . . . . . . . . . . . . . . . . 86 v 3.5.2.1 Effect of particle spacing . . . . . . . . . . . . . . . . . 86 3.5.2.2 Maximum concentration . . . . . . . . . . . . . . . . . 89 3.5.2.3 Negative concentrations . . . . . . . . . . . . . . . . . 89 3.5.2.4 Effect of ratio between smoothing length and particle spacing . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.5.2.5 Effect of anisotropy ratio . . . . . . . . . . . . . . . . . 93 3.5.2.6 Effect of kernel function . . . . . . . . . . . . . . . . . 93 3.5.2.7 Effect of velocity orientation . . . . . . . . . . . . . . . 96 3.5.3 Irregularly Spaced Particles . . . . . . . . . . . . . . . . . . . . 97 3.5.3.1 Isotropic case . . . . . . . . . . . . . . . . . . . . . . . 99 3.5.3.2 Anisotropic case . . . . . . . . . . . . . . . . . . . . . 104 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 3.7 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4 A Multidimensional Streamline-Based Method to Simulate Reactive Solute Transport in Heterogeneous Porous Media . . . . . . . . . . . 111 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.2 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.2.1 Governing Equation . . . . . . . . . . . . . . . . . . . . . . . . 116 4.2.2 Streamline Formulation . . . . . . . . . . . . . . . . . . . . . . . 117 4.2.3 Numerical Approximation . . . . . . . . . . . . . . . . . . . . . 118 4.2.3.1 Advection along streamlines . . . . . . . . . . . . . . . 119 4.2.3.2 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.3 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.3.1 Streamline Tracing . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.3.2 Time Integration . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.3.3 Advection Solution . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.3.4 MC-SPH Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 125 4.3.4.1 SPH kernel . . . . . . . . . . . . . . . . . . . . . . . . 125 4.3.4.2 Neighbor search . . . . . . . . . . . . . . . . . . . . . . 125 4.3.4.3 Time integration . . . . . . . . . . . . . . . . . . . . . 126 4.3.5 Longitudinal Dispersion . . . . . . . . . . . . . . . . . . . . . . 127 4.3.5.1 Interface coefficients . . . . . . . . . . . . . . . . . . . 127 4.3.5.2 Time integration . . . . . . . . . . . . . . . . . . . . . 127 4.4 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 vi 4.4.1 Example 1: Continuous Solute Release in Uniform Flow . . . . . 130 4.4.2 Example 2: Quarter Five-Spot in Heterogeneous Medium . . . . 132 4.4.2.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 4.4.2.2 Simulated concentrations . . . . . . . . . . . . . . . . 134 4.4.2.3 Numerical oscillations . . . . . . . . . . . . . . . . . . 140 4.4.2.4 Performance comparison . . . . . . . . . . . . . . . . . 142 4.4.3 Example 3: Quarter Five-Spot in Heterogeneous Medium with Rate-Limited Sorption . . . . . . . . . . . . . . . . . . . . . . . 143 4.4.4 Example 4: Natural Biodegradation in Three-dimensional Hetero- geneous Porous Media . . . . . . . . . . . . . . . . . . . . . . . 146 4.4.4.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 4.4.4.2 Simulated concentrations . . . . . . . . . . . . . . . . 149 4.4.4.3 Breakthrough curves . . . . . . . . . . . . . . . . . . . 156 4.4.4.4 Performance . . . . . . . . . . . . . . . . . . . . . . . . 159 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 4.6 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 5.1 Limitations of Proposed Numerical Schemes . . . . . . . . . . . . . . . 173 5.2 General Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 5.3 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 5.4 Final Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 5.5 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 Appendix A Derivation of SPH Approximation for Isotropic Dispersion 181 Appendix B Derivation of SPH Approximation for Second Order Deriva- tives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 Appendix C Random Walk Particle Method . . . . . . . . . . . . . . . 188 Appendix D Streamline Tracing . . . . . . . . . . . . . . . . . . . . . . . 191 vii List of Tables 1.1 Estimated dispersivity values from field- and laboratory-scale experiments. 9 2.1 Comparison of MC-SPH and RWPT methods. . . . . . . . . . . . . . 44 2.2 Parameters used in RWPT simulations. . . . . . . . . . . . . . . . . . . 45 2.3 Parameters used in SPH simulations. . . . . . . . . . . . . . . . . . . . 45 2.4 Parameter and results of flow model. . . . . . . . . . . . . . . . . . . . 54 2.5 Parameter values used in transport model. . . . . . . . . . . . . . . . . 55 3.1 Parameters used in all simulations. . . . . . . . . . . . . . . . . . . . . 82 3.2 Parameters used to define different simulation scenarios to evaluate ap- proximations for anisotropic dispersion. . . . . . . . . . . . . . . . . . 85 3.3 Definition of different runs used to study convergence properties. . . . . 85 3.4 Normalized error for different SPH kernels and PSE cutoff functions. . 96 3.5 Error versus flow velocity direction. . . . . . . . . . . . . . . . . . . . . 97 4.1 Dispersivity and equivalent longitudinal (PeL) and transverse (PeT ) grid Péclet values used in Example 1. . . . . . . . . . . . . . . . . . . . . . 131 4.2 Parameters used in MOC simulations. . . . . . . . . . . . . . . . . . . . 133 4.3 Number of nodes or cells, time step size and number of time steps used in simulations of Example 2. . . . . . . . . . . . . . . . . . . . . . . . . . 134 4.4 Dispersivity and equivalent longitudinal (PeL) and transverse (PeT ) grid Péclet values used in Example 2. . . . . . . . . . . . . . . . . . . . . . 134 4.5 Normalized minimum simulated concentration values for Example 2. . . 142 4.6 Normalized maximum simulated concentration values for Example 2. . 142 4.7 Normalized CPU time required to simulate Example 2 for different sce- narios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 4.8 Parameters of the rate-limited sorption model used in Example 3. . . . 144 4.9 Definition of three scenarios simulated in Example 4. . . . . . . . . . . 147 4.10 Spatial and temporal discretizations used to simulate Example 4. . . . 149 viii 4.11 Normalized CPU time required to simulate Example 4 for the two scenar- ios than include biodegradation. . . . . . . . . . . . . . . . . . . . . . . 162 D.1 Comparison of the performance of Pollock’s and explicit adaptive algo- rithm to trace streamlines in heterogeneous quarter five-spot problem. . 204 D.2 Comparison of alternative seed distributions to trace streamlines in the heterogeneous quarter five-spot problem. . . . . . . . . . . . . . . . . 209 ix List of Figures 1.1 Groundwater pollution due to tailings infiltration. . . . . . . . . . . . 2 1.2 Vertical cross-section of a synthetically generated aquifer using the esti- mated statistics of the sandy aquifer at Canadian Forces Base, Borden, Ontario. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Solute plume migration in vertical cross-section. . . . . . . . . . . . . 8 1.4 Effect of advection and local-scale dispersion on solute concentration. 11 1.5 Subgrid-scale segregation and cell averaged concentration values. . . . 14 2.1 SPH kernels and derivatives. . . . . . . . . . . . . . . . . . . . . . . . 34 2.2 Error versus smoothing length. . . . . . . . . . . . . . . . . . . . . . . 43 2.3 Initial distribution of particles and solute concentration in RWPT sim- ulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.4 Initial solute concentration distribution in MC-SPH simulations. . . . 47 2.5 Normalized CPU time versus number of particles in RWPT simulations. 49 2.6 Normalized CPU time versus total number of particles and average num- ber of particles per kernel support volume in MC-SOH simulations. . 50 2.7 Maximum concentration versus time in MC-SPH and RWPT simulations. 52 2.8 Normalized global error versus total number of particles and normalized CPU time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 2.9 Domain dimensions, square initial plume, and breakthrough observation points P1 and P2 along centerline. . . . . . . . . . . . . . . . . . . . . 55 2.10 Spatial concentration distribution for TVD, HMOC and SPH simula- tions at τ = Ut/IY = 62 for Pe =∞. . . . . . . . . . . . . . . . . . . 58 2.11 Concentration versus accumulated distance along centerline at dimen- sionless time τ = Ut/IY = 62. . . . . . . . . . . . . . . . . . . . . . . 59 2.12 Breakthrough curve at point P1 located 26IY downstream from initial plume center. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 x 2.13 Breakthrough curve at point P2 located 42IY downstream from initial plume center. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.14 Dimensionless maximum concentration versus dimensionless time. . . 63 2.15 Mean concentration versus dimensionless time. . . . . . . . . . . . . . 65 3.1 Sum of dispersion components as function of flow velocity direction. . 79 3.2 SPH kernels and PSE cutoff functions. . . . . . . . . . . . . . . . . . 84 3.3 Error E2 as function of particle or grid spacing for equispaced particles. 87 3.4 Normalized error E∞ as function of particle or grid spacing for equis- paced particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.5 Difference between maximum concentration values of numerical and an- alytical solutions as function of time for equispaced particles. . . . . . 90 3.6 Concentration distribution for equispaced particles. . . . . . . . . . . 91 3.7 Difference between analytical and numerical solutions for equispaced particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.8 Error E2 versus the ratio between smoothing length or core size and particle spacing, γ = h/∆x. . . . . . . . . . . . . . . . . . . . . . . . 94 3.9 Error as function of the anisotropy ratio αT/αL. . . . . . . . . . . . . 95 3.10 Particle distortion due to flow velocity. . . . . . . . . . . . . . . . . . 98 3.11 Particle locations considering random and quasi-random distributions. 99 3.12 E2 error versus average particle spacing using equispaced, random, and quasi-random particle distributions. . . . . . . . . . . . . . . . . . . . 100 3.13 Normalized E∞ error versus average particle spacing using equispaced, random, and quasi-random particle distributions. . . . . . . . . . . . 102 3.14 Difference between maximum concentration values of analytical and nu- merical solutions as function of time for random and quasi-random par- ticle distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.15 Error as function of average particle spacing for αT/αL = 0.01 using randomly and quasi-randomly distributed particles. . . . . . . . . . . 105 3.16 Concentration distribution after 300 time steps for αT/αL = 0.01 and quasi-randomly distributed particles. . . . . . . . . . . . . . . . . . . 106 4.1 Meshless MC-SPH method. . . . . . . . . . . . . . . . . . . . . . . . . 113 4.2 Hybrid streamline-SPH method. . . . . . . . . . . . . . . . . . . . . . 115 xi 4.3 Flow oriented coordinate system. . . . . . . . . . . . . . . . . . . . . 118 4.4 Overall solution approach implemented in streamline-based simulator. 129 4.5 Comparison of simulated concentrations for Example 1. . . . . . . . . 131 4.6 Spatial distribution of the natural logarithm of the hydraulic conductiv- ity and streamlines in Example 2. . . . . . . . . . . . . . . . . . . . . 132 4.7 Simulated concentration values after injection of 0.4 pore volume of con- taminated fluid for Example 2. . . . . . . . . . . . . . . . . . . . . . . 135 4.8 Breakthrough curves at observation point P1 in Example 4. . . . . . . 137 4.9 Breakthrough curves at observation point P2 in Example 4. . . . . . . 138 4.10 Comparison of simulated breakthrough curves at observation point P1 for Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 4.11 Comparison of simulated breakthrough curves at observation point P2 for Example 2, (a) streamline simulator and (b) MOC solver using fine grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 4.12 Numerical oscillations in simulated concentrations for Example 2. . . 141 4.13 Breakthrough at observation point P2 for different mass transfer (β) and partition (Kd) coefficients considered in Example 3 . . . . . . . . . . 145 4.14 Spatial distribution of natural logarithm of hydraulic conductivity and flow velocity magnitude used in Example 4 . . . . . . . . . . . . . . . 148 4.15 Simulated concentrations at nodes along streamlines after 10,000 days since the initial release of BTEX for the scenario that includes advective transport with biodegradation in Example 4. . . . . . . . . . . . . . . 151 4.16 Simulated concentrations at nodes along streamlines after 10,000 days since the initial release for the scenario that includes advection, disper- sion and biodegradation in Example 4. . . . . . . . . . . . . . . . . . 152 4.17 imulated concentration values for Example 4 at vertical plane defined by y=11.5 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 4.18 Simulated concentration values for Example 4 at horizontal plane defined by z=2.5 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 4.19 Comparison of simulated concentration values for Example 4 at horizon- tal plane defined by z=2.5 m. . . . . . . . . . . . . . . . . . . . . . . 155 4.20 Simulated concentration values for Example 4 along the profile parallel to y direction at coordinates x=35 m and z=2.5 m. . . . . . . . . . . 157 xii 4.21 Location of contaminant source and the two observation wells in Exam- ple 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 4.22 Simulated BTEX concentration versus time for Example 4 assuming advective transport only. . . . . . . . . . . . . . . . . . . . . . . . . . 160 4.23 Simulated BTEX concentration versus time for Example 4 assuming advective transport and biodegradation. . . . . . . . . . . . . . . . . . 161 4.24 Distribution of cells according to the flow velocity magnitude. . . . . 163 4.25 Distribution of (a) nodes along streamlines according to the flow veloc- ity magnitude and (b) number of streamlines based on the maximum velocity magnitude along individual streamlines. . . . . . . . . . . . . 164 4.26 Spatial distribution of nodes and number of nodes according to the num- ber of neighboring nodes that contribute to the SPH summation to ap- proximate dispersion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 D.1 Scehamatic of Pollock’s particle tracking method. . . . . . . . . . . . 193 D.2 Schematic of explicit particle tracking method. . . . . . . . . . . . . . 198 D.3 Comparison of Pollock’s and explicit integration methods for the homo- geneous quarter five-spot problem. . . . . . . . . . . . . . . . . . . . 200 D.4 Comparison of time of flight and arc length computed with Pollock’s and explicit integration methods for homogeneous quarter five-spot problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201 D.5 Spatial distribution of hydraulic conductivity used in heterogeneous quater five-spot problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201 D.6 Comparison of Pollock’s and explicit integration methods for the het- erogeneous quarter five-spot problem. . . . . . . . . . . . . . . . . . 202 D.7 Comparison of time of flight and arc length computed with Pollock’s and explicit integration methods for heterogeneous quarter five-spot problem. 203 D.8 Streamline distribution in heterogeneous quater five-spot problem. . . 207 D.9 Comparison of two alternative streamline distributions for the homoge- neous quater five-spot problem. . . . . . . . . . . . . . . . . . . . . . 207 D.10 Comparison of four alternative streamline distributions for the hetero- geneous quater five-spot problem. . . . . . . . . . . . . . . . . . . . . 208 D.11 Schematic of flow problem considering low permeability inclusion. . . 210 xiii D.12 Comparison of Pollock’s and time of flight streamline discretization ap- proaches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 D.13 Comparison of Pollock’s and arc length streamline discretization ap- proaches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212 xiv Acknowledgments I especially want to thank my advisor, Prof. Roger Beckie, who gave me freedom and constant encouragement to pursue my research interests. He has been an excellent mentor not only on the technical aspects of this research, but also on the many other aspects of the academic world. I feel fortunate I have had the opportunity to collaborate with such great person. I also appreciate the interesting conversations I had with the others members of my thesis committee, Prof. Leslie Smith and Prof. Ulrich Mayer. Prof. Albert Valocchi and Dr. Marco Massabó also provided comments that improved the text and help to clarify some of the points presented in this thesis. I am grateful to all the people that directly or indirectly have funded my education and research through the years. My PhD education and research have been founded through a University of British Columbia Graduate Fellowship and a NSERC Discovery Grant awarded to Prof. R. Beckie. I am also grateful to have received the Thomas and Marguerite MacKay Memorial Scholarship and the Hugh Nasmith Graduate Scholarship while I was a graduate student at UBC. I am deeply indebted with thousands of programmers who have generously donated their work to create free software. This thesis would not have been possible without such great tools. Finally, I wish to thank my family who has been a constant source of support and encouragement during these long years. To Françoise, Clara and Matias, I will never forget all the sacrifices you have made to allow me to pursue my dreams. xv To Françoise, Clara and Matias xvi Co-authorship Statement This thesis has been prepared as a collection of manuscripts, either published, submitted or in preparation, that have been co-authored with individuals other than myself. In each case I am the first author and have conducted all the research work and manuscript preparation. The specific objectives of each chapter and research approach are based on my initiative in consultation with my thesis supervisor Prof. Roger Beckie. The research work, including derivation of numerical approximations, code development and numerical simulations, with the exception of the one-dimensional simulations presented in Chapter 2, have been entirely done by myself with support from the co-authors as outlined below. Dr. Marco Massabó from CIMA Research Foundation in Savona, Italy; provided insight- ful comments and collaboration on the interpretation of the smoothed particle hydrody- namics (SPH) method in the context of simulations of subsurface solute transport. In addition, Dr. Massabó provided the results of the one-dimensional SPH simulations in- cluded in Chapter 2. Prof. Beckie, Dr. Massabó, and Prof. Valocchi from the University of Illinois at Urbana-Champaign, USA, also assisted with corrections and suggestions to improve the editing of the chapters which they were involved with. xvii Chapter 1 Introduction 1.1 Motivation Groundwater pollution has become a serious problem during recent decades. The release and infiltration into aquifers of pesticides, organic compounds, contaminants of biolog- ical origin and nuclear waste, among others, pose a risk for the human health and the environment. In 1996, the U.S. Environmental Protection Agency (EPA) identified 217,000 contam- inated sites in the U.S. that required mitigation action. Another 300,000 sites were reported cleaned up or were found to no require mitigation. It was estimated that sev- eral thousands of those sites were polluted with highly radioactive nuclear waste and would require coordinated mitigation actions for several decades before been declared cleaned up. In addition, it was estimated that there were between 130,000 to 450,000 additional sites that could potentially require some mitigation action (EPA, 1996). In Canada, the Federal Contaminated Site Inventory (FCS, retrieved on August 10, 2009) lists 3,208 sites that have at least one substance in the groundwater that occurs at concentrations above natural levels and that pose an immediate or long-term hazard to human health or the environment. Those sites includes only the small proportion of cases for which the Government of Canada has accepted some or all financial responsibility. Figure 1.1 shows a common example of groundwater contamination due to tailings infil- tration from a tailings impoundment. Once the tailings plume reaches the water table it migrates carried by the regional groundwater flow and it can eventually impact the 1 Figure 1.1: Groundwater pollution due to tailings infiltration. Once the tailings plume reaches the water table it can migrate up to several kilometers downstream from the contaminant source carried by the regional ground- water flow. water quality of wells located up to several kilometers downstream from the contaminant source. This example demonstrates the potential large spatial scale of problems related to the pollution of natural aquifers. Because of the scale of the problem, it is essential to find the most effective and financially sound mitigation actions. Possible mitigation actions include: natural attenuation (ei- ther due to dilution or natural biodegradation), passive containment, and active cleanup measures such as pump and treat and enhanced bioremediation (EPA, 1996). The selec- tion of the most effective action requires a good understanding of the physical, chemical and biological processes that govern the migration and transformation of contaminants in the subsurface. At the same time, scientists and engineers who are involved in the remediation of contaminated sites are interested in finding answers to questions like: • How long will the contaminant plume take to reach a well or a river? • What will be the contaminant concentration at a given location and time? • Will dilution due to the advection and dispersion of the contaminant plume be enough to decrease the contaminant concentration to acceptable levels within a reasonable time frame? 2 • Will the contaminant of concern be retarded with respect to the groundwater flow? • Will biodegradation be an effective process to transform and remove the contami- nant from the groundwater? Reactive transport modeling has emerged during the last decades as an instrument to answer these practical questions and as a tool to integrate fundamental knowledge of the complex processes that control flow, transport and chemical reactions in porous media (Steefel et al., 2005). In this thesis, new methods for the simulation of conservative and reactive transport are investigated. The methods focus upon accurate and efficient solution of physical transport processes in heterogeneous porous media. This research addresses problems in existing formulations that lead to, for example, negative concentrations, which are particularly problematic because of non-linear chemical reaction rates that are common in reactive transport simulations. 1.2 Reactive Transport Modeling 1.2.1 Definition Reactive transport in natural porous media is a broad term that is used to refer to complex physical and chemical processes that occur at disparate spatial and temporal scales and that involve fluid flow, mass transport and chemical reactions in the subsurface (Steefel et al., 2005). The interaction between transport, and reactions is complex. On one hand, reactions such as mineral precipitation and dissolution can change the porosity and permeability of a porous medium, hence, affect the fluid flow and transport properties of the medium. On the other hand, mass transport plays a key role in enabling reactions because it provides the driving force to perturb a chemical system out of equilibrium by transport- ing and mixing reactants and, because it sets a characteristic time scale during which reactions can take place (Valocchi, 1985). The most important transport processes for enabling reactions are advection, molecular diffusion, and mechanical dispersion (Steefel and MacQuarrie, 1996; De Simoni et al., 2005; Steefel and Maher , 2009). 3 1.2.2 Conceptual Model The migration and transformation of contaminants in the subsurface is the result of fluid and transport processes that occur at length scales of a single or few pores and of chemical reactions that happen at even smaller scales (molecular scale). Although, pore-scale modeling has gained considerable attraction during recent years, their use is limited to very small scale problems where they have provided insightful understanding of basic mechanisms. However, pore-scale modeling is not appropriate for studying practical field-scale problems that occur at much larger scales (Steefel et al., 2005). Therefore, most current conservative and reactive transport models used to simulate field-scale problems are based on a continuum representation of porous media, such that the system properties are averaged over a representative elementary volume (REV) with length scale equal to many pore lengths (Bear , 1988; Steefel et al., 2005). Thus, the REV scale, often called Darcy’s scale or local-scale, defines the spatial scale or volume size where fluid velocity, transport properties, concentrations and reactions rates are computed. In what follows we will concentrate our discussion on the use of reactive transport model- ing in hydrogeology. Furthermore, we will assume isothermal saturated groundwater flow with constant density and negligible effect of reactions on flow and transport properties. 1.2.3 Mathematical Model The groundwater specific discharge, q, can be calculated using Darcy’s law, q = −K∇φ (1.1) where K [L/T] is the hydraulic conductivity tensor and φ [L] is the hydraulic head. The two main transport mechanisms at the REV scale are: advection, which involves the movement of the solute with the flow; and hydrodynamic or local-scale dispersion that includes molecular diffusion and mechanical dispersion due to variations of the flow velocity at the pore-scale (Bear , 1988). Thus, reactive solute transport at the local-scale is modeled by a system of partial differential equations, which is given as follows for the case of constant porosity (Bear , 1988; Steefel and MacQuarrie, 1996): 4 ∂Ck ∂t = ∇ · (D∇Ck)−∇ · (vCk) +Rk(c) k = 1, . . . ,m ∂Ck ∂t = Rk(c) k = m+ 1, . . .M (1.2) where Ck [M/L3] is the solute concentration of species or component k, D [L2/T] is the hydrodynamic dispersion tensor, v = q/η [L/T] is the pore water velocity, η is the porosity of the medium, Rk [M/L3/T] is the total reaction rate for species or component k, c = (C1, . . . , CM) is the concentration vector,m is the number of species or components in the aqueous (mobile) phase, and M is the total number of species or components. The most common expression to compute the coefficients of the dispersion tensor D for an isotropic porous medium considering a Cartesian coordinate system is (Bear , 1988) Dij = (αT q +Dm) δij + (αL − αT ) vivj q (1.3) where Dm [L2/T] is the molecular diffusivity, δij is Kronecker’s delta, q = |v| [L/T] is the magnitude of the pore water velocity, and αL and αT [L] are the longitudinal and transverse dispersivity of the medium, respectively. Alternative expressions for the dispersion tensor components include different transverse dispersivities for the horizontal and vertical directions (Burnett and Frind, 1987; Lichtner et al., 2002), however, those models are less commonly used. The reaction term in (1.2) may include homogeneous reactions that occur in a single phase or heterogeneous reactions that include constituents in more than one phase, e.g. sorption which includes the solid and aqueous phases (Rubin, 1983; Mayer et al., 2002). From a practical point of view, sorption and biodegradation are the two most relevant reactions in groundwater, because of their role in retarding the migration of heavy metals and in the transformation of hydrocarbons, respectively; which are the two most common substances in contaminated aquifers (EPA, 1996). 5 1.3 Numerical Solution The system of equations in (1.2) corresponds to a set of non-linear partial differential equations (PDE) and, in general, must be integrated numerically. There are two different numerical approaches to integrate the system of PDEs: a fully-implicit approximation or an operator splitting approach (Yeh and Tripathi, 1989; Steefel and MacQuarrie, 1996). The simplest one is based on an operator splitting formulation that allows decoupling of the transport and reaction terms. Then, equation (1.2) is split into two terms rep- resenting transport and chemical reactions. Usually, the transport component, which corresponds to a linear PDE, is solved first, then the concentrations computed as result of the transport step are used as initial conditions to compute the solution of the non- linear set of ordinary differential equations that represent chemical reactions. Variations of the operator splitting approach include schemes that iterate between the solutions of the transport and reaction terms or that switch the order of the evaluation (Steefel and MacQuarrie, 1996). In the discussion that follows we will assume that an operator splitting approach is used to evaluate (1.2) and we will focus our analysis on the numerical solution of the transport term. However, we must emphasize that the analysis presented below is also valid for fully implicit implementations. 1.3.1 Particularities of Flow and Transport in Porous Media The numerical integration of the advection-dispersion-reaction (ADR) equation, (1.2), presents some unique challenges. First, because of the geological origin of aquifers, hydraulic conductivity varies by sev- eral orders of magnitudes within relatively short distances. For example, Sudicky (1986) reports differences of more than thirty times between hydraulic conductivity values sepa- rated by few centimeters in a shallow sandy aquifer at Canadian Forces Base in Borden, Ontario, Canada. Similar variations were also observed at the Macrodispersion Experi- ment (MADE) site in Mississippi, USA, where local-scale hydraulic conductivity values varied by more than four orders of magnitude within an area of approximately 250 x 300 m (Boggs et al., 1992; Zheng and Gorelick, 2003). Large variations in hydraulic conductivity produce not only important variations in flow velocity magnitude, but also in the direction of the flow (Sudicky, 1986). For example, 6 Figure 1.2: Vertical cross-section of a synthetically generated aquifer using the esti- mated statistics of the sandy aquifer at Canadian Forces Base, Borden, Ontario (Mackay et al., 1986; Sudicky, 1986; Freyberg, 1986). Natural logarithm of the hydraulic conductivity (top) and simulated groundwater velocity (bottom). The cross-section is 32 m long and 8 m high. Large variations of hydraulic conductivity within short distances result in im- portant variation of the magnitude and direction (arrows) of the flow velocity. Figure 1.2 shows a vertical cross-section of a synthetically generated aquifer using the estimated statistics of the sandy aquifer at the Borden site (Mackay et al., 1986; Sudicky, 1986; Freyberg, 1986). The figure also shows simulated groundwater velocities, that demonstrate that variations of hydraulic conductivity within short distances result in important variations of the magnitude and direction of the flow velocity. Because of the heterogeneity of the flow velocity, adjacent fluid parcels may travel at very different velocities. Such variations in travel time result in stretching and spreading of the contaminant plume, which manifests as large variations of concentration within short distances even at the order of few centimeters (Mackay et al., 1986;Molz and Widdowson, 7 1988; Benson and Meerschaert, 2008). Figure 1.3 shows how variations in velocity can produce dramatic changes in the shape of a initially regular plume within short travel distances. Figure 1.3: Solute plume migration in vertical cross-section shown in Figure 1.2. An idealized initial rectangular solute plume (red rectangle on top figure) migrates carried by the flow velocity. Because of the heterogeneity of the flow field, the shape of the plume becomes very irregular after the center of mass of the plume has travelled approximately 10 m (bottom). Second, local-scale dispersion is anisotropic and it is much more important in the direction of the flow than in the transverse directions. Table 1.1 lists estimated dispersivity values based on data collected in field- and laboratory-scale experiments. While longitudinal dispersivity ranges between 3 · 10−2 and 5 · 10−1 m, typical transverse dispersivity is of the order of 10−3 m, which is equivalent to transverse dispersion coefficients of the order of molecular diffusion for typical pore water velocities (Cirpka et al., 2006). Note that experimental measurements of local-scale dispersion are difficult, thus most estimated dispersivity values are based on spatial moment analysis, which measures the spreading of the plume, or effective parameters that measure the combined effect of flow heterogeneity and local-scale dispersion (Jose et al., 2004). Therefore, most of the values presented in 8 Table 1.1 represent upper bounds for local-scale longitudinal and transverse dispersivities. αL (m) αT (m) Comments 3 · 10−2 − 5 · 10−1 5 · 10−4 − 1 · 10−3 Based on data collected during field-scale tracer test in a shallow unconfined sand and gravel aquifer on Cape Cod, Massachusetts, USA (Hess et al., 2002) 4 · 10−1 4 · 10−2 Based on spatial moment analysis of data collected during the large-scale experiment conducted at the Borden site in Ontario, Canada (Freyberg, 1986) – 1 · 10−6 − 2 · 10−4 Calculated from reactive plume lengths in laboratory-scale experiments with homogeneous materials (Cirpka et al., 2006) 4 · 10−1 2 · 10−3 Effective dispersivities obtained from analysis of breakthrough curves in a 14 m long sandbox filled with four different types of silica sand (Jose et al., 2004) Table 1.1: Estimated dispersivity values from field- and laboratory-scale experiments. The heterogeneity of the flow velocity field together with the small magnitude of local- scale dispersion coefficients in natural porous media, have important consequences for key transport processes such as spreading, dilution and mixing. Kitanidis (1994) discusses the difference between the spreading and dilution of a solute plume. Spreading is defined as the stretching of the plume and can be measured as the rate of change of the second central spatial moment. Dilution is the process by which the initial solute mass is distributed in an increasing volume. Mixing is the result of the combined action of the stretching and folding of material lines of the plume, and the mass exchange due to local-scale dispersion (Weeks and Sposito, 1998). While spreading is produced by the spatial velocity variability, dilution and mixing are due to the combined action of the heterogeneity of the flow field and most importantly local-scale dispersion. Because of the relatively small magnitude of the local-scale dispersion and in absence of sorption, mixing between a contaminant and other chemical species present in the 9 natural groundwater occurs in a narrow zone located along the irregular edges of the contaminant plume (Oya and Valocchi, 1998). Figure 1.4 shows the difference between spreading and mixing using results of simulations included in Section 2.3.3. The figures show simulated concentrations for two conservative transport scenarios with different Péclet number (Pe): advection-only (Pe =∞) and for advection and dispersion (Pe = 200). Under a purely-advective scenario the solute plume becomes irregular due to variations in velocity, however initial concentration values do not change and there is a sharp interface between fluid zones with and without solute. When local-dispersion is included, there is solute mass transfer in areas of high concentrations near the plume edges that results in a thin area, relative to the typical length of the local-scale heterogeneity, with concentrations lower than the initial value. In many field situations the zone with lower concentrations would correspond to a mixing area where the solute and the natural groundwater mix enabling chemical reactions. For example, natural and enhanced attenuation of organic contaminants occur in a narrow zone close to the plume boundaries where the contaminant (substrate) and the electron acceptor (e.g. oxygen) mix (Oya and Valocchi, 1998; Cirpka et al., 1999b; Ham et al., 2004). In the next section, we will argue that the characteristics of the flow and solute transport process described above must be considered in the selection of numerical methods to simulate conservative and reactive solute transport in groundwater. 1.3.2 Numerical Methods We start this section by reviewing some of the most common numerical schemes that are used to simulate conservative and reactive solute transport. 1.3.2.1 Mesh-based numerical methods This category includes finite difference, finite volume or finite element methods. In finite difference and finite volume schemes each grid cell defines a new control volume within which parameters and variables are considered constant (Steefel and MacQuarrie, 1996). Mesh-based methods are relatively easy to implement, have convergence, stability and accuracy properties that are well understood, and it is possible to develop formulations that are mass conservative. Because of those characteristics, mesh-based approximations 10 Figure 1.4: Effect of advection and local-scale dispersion on solute concentration. Fig- ure shows a small part of the simulated solute plume in Section 2.3.3 for Péclet number, Pe = ∞ (top) and Pe = 200 (bottom). Advection only affects the shape of a contaminant plume (top). If local-dispersion is in- cluded, mass transfer occurs in areas of high concentration gradients and a mixing zone develops around the plume edges (bottom). The mixing zone is critical to enable some chemical reactions such as biodegradation. 11 are used in most reactive transport packages (e.g. Pruess, 1991;White et al., 1995;Mayer et al., 2002; Mills et al., 2007). Low-order mesh-based schemes to approximate advection, e.g. upstream finite differ- ence, introduce large amount of artificial diffusion and mixing (Steefel and MacQuarrie, 1996). Therefore, alternative schemes based on high-order approximations that were first developed to simulate problems in computational fluid dynamics, have been adopted by reactive transport modelers. Examples of high-order schemes are flux-corrected transport (FCT) methods (Boris and Book, 1973; Zalesak, 1979), which combine high- and low- order schemes; and total variation diminishing (TVD) schemes (Harten and Lax , 1984; Yee et al., 1985; Cox and Nishikawa, 1991). 1.3.2.2 Hybrid Eulerian-Lagrangian methods Hybrid schemes are based on the same general concept, the use of particles to handle advection and a grid-based method to handle dispersion. Each time step is split into two sub-steps. First, changes in concentrations due to advection are computed by forward or backward particle tracking. Then, concentration values are interpolated onto a grid. Next, grid concentration values are used to solve for dispersion, and eventually reactions, using some traditional mesh-based solver. Multiple variations of this approach exist depending on the interpolation methods and tracking algorithm. Examples of these kinds of methods are: hybrid Eulerian-Lagrangian methods (Neuman, 1981, 1984), method of characteristics (MOC) and hybrid method of characteristics (HMOC) in the MT3DMS package (Zheng and Wang, 1999) and MOC3D (Konikow et al., 1996), and Eulerian- Lagrangian localized adjoint methods (ELLAM) (Celia et al., 1990; Russell and Celia, 2002). 1.3.2.3 Random walk particle tracking methods Random walk particle-tracking methods (RWPT) have long been used to simulate con- servative solute transport in porous media (Ahlstrom et al., 1977; Pickens and Grisak, 1981; Tompson and Gelhar , 1990; Tompson, 1993). In this type of model, solute mass is distributed among a set of particles that move carried by the flow velocity and by a random drift that models dispersive transport (Delay et al., 2005; Salamon et al., 2006). Solute concentrations are estimated by averaging the mass contained in the particles 12 found in some specified volume. Therefore, concentration values depend upon the total number of particles, size of the averaging volume, and spatial particle distribution. The popularity of RWPT methods is due to its natural capacity to accurately simulate advection and ease of implementation. Because of its advantages RWPT has become the de facto standard method in numerical studies of plume spreading and dilution (e.g. see Delay et al., 2005; Salamon et al., 2006, and references therein). However, RWPT methods are less attractive for the simulation of reactive transport because: (i) it is difficult to simulate general heterogeneous reactions that include the solid phase, (ii) a very large number of particles is required to obtain an accurate estimation of low concentration values that is crucial to approximate reactions that occur in the mixing zone along the plume edges, and (iii) simulations that include multiple species need a large number of particles to track individual species. 1.3.3 Limitations of Current Numerical Methods According to Steefel and MacQuarrie (1996) there are three main properties that a nu- merical method must satisfy to be used in reactive transport simulations: (i) accuracy in space and time, which includes minimizing numerical diffusion and mass conserva- tion errors, (ii) monotonicity, which means avoiding spurious oscillations (e.g., negative concentrations); and (iii) computational efficiency. Next, we evaluate current numerical methods based on those three criteria. 1.3.3.1 Accuracy Mesh-based numerical methods, including high resolution methods, have problems to ac- curately simulate multidimensional advection-dominated transport because of numerical dispersion that results in excessive artificial mixing, dilution, and overestimation of reac- tions rates (Steefel and MacQuarrie, 1996; Cirpka et al., 1999; Zheng and Wang, 1999). Numerical dispersion is more important when the grid is not aligned with the direction of the flow (Frind et al., 1987), which is always the case in non-uniform flows as found in heterogeneous porous media. On the other hand, hybrid schemes that require interpolating concentrations to a back- ground grid also introduce numerical dispersion even if they provide a very accurate solution for advection. In the case of RWPT methods, concentration values can only 13 A! B! A + B! Figure 1.5: Subgrid-scale segregation and cell averaged concentration values. Species A and B are physically segregated at the subgrid-scale (left). However, they would appear well-mixed at the cell scale (right). If A and B are two reactants in a chemical reaction, then numerical simulations based on cell averaged concentrations would overestimate the reaction rate. be obtained after averaging the mass of particles over some control volume, which also results in numerical mixing. The use of cell average concentration values can introduce large errors in the estimation of dilution and reaction rates. For example, Figure 1.5 shows two species, A and B, which are physically segregated, however, they would appear well-mixed in numerical simulations that use cell averaged concentration values. If A and B are reactants in a chemical reaction, then the simulated reaction rate would overestimate the real reaction rate (which is zero, in this case). The errors due to numerical dispersion and cell averaging are smaller for larger values of local-scale dispersion, because concentration values within the cell volume are smoothed out by dispersion (Steefel and MacQuarrie, 1996). However, as discussed earlier, local- scale dispersion in porous media, particularly in the transverse direction, is small and its effect to smooth out concentration fluctuations is limited. Therefore, sub-grid scale concentration fluctuations are important and the use of cell averaged concentration val- ues is an important source of error in conservative and reactive transport simulations (Frind and Germain, 1986; Frind et al., 1987; Molz and Widdowson, 1988; Benson and Meerschaert, 2008). In evaluating the accuracy of a numerical method for simulating reactive transport, it 14 is also important to keep in mind that any error in the solution of the transport com- ponent can be greatly amplified by non-linear reactions. Thus, methods that perform acceptably well to simulate conservative solute transport, can produce large errors when chemical reactions are included (Steefel and MacQuarrie, 1996). For example, Cirpka et al. (1999b) demonstrated that small amounts of numerical dispersion in simulations of biodegradation controlled by transverse mixing simulated using a high-order finite volume method, can result in larger errors in the estimation of reaction rates and contaminant mass removal. 1.3.3.2 Monotonicity Spurious oscillations in simulations of conservative solute transport arise due to the use of non-linear high-order methods to control numerical dispersion (Steefel and MacQuar- rie, 1996; Cirpka et al., 1999) and numerical approximations of the off-diagonal entries (“cross-terms”) in the local-dispersion tensor (Herrera and Valocchi, 2006). Multidimen- sional high-order mesh-based solvers for advection based on the FCT and TVD schemes, which are supposed to suppress numerical oscillations, often result in small oscillations (Steefel and MacQuarrie, 1996; Herrera and Valocchi, 2006). Numerical oscillations that arise from the solution of parabolic or elliptic PDEs that include mixed derivatives or “cross-terms” are a well known problem and have been the subject of many research efforts in recent years (e.g. Nordbotten and Aavatsmark, 2005; Le Potier , 2005b; Mlacnik and Durlofsky, 2006; Edwards and Zheng, 2008; Yuan and Sheng, 2008; Lipnikov et al., 2009). To this day, no single solution provides a scheme that can be used in general scenarios. Although, small oscillations can be usually tolerated in conservative transport simula- tions, they are unacceptable in reactive transport simulations because they can result in unstable solutions in presence of non-linear chemistry. For example, Steefel and Mac- Quarrie (1996) discusses the effect of small oscillations in a problem involving organic car- bon degradation via sulfate reduction coupled to two equilibrium dissolution-precipitation reactions. They simulate solute transport using a high-order FCT method that intro- duces spurious oscillations that do not produce problems in a tracer simulation, but that produce unstable results when chemical reactions are included. 15 1.3.3.3 Performance Numerical simulations of reactive transport in porous media are computationally de- manding. Although, the increasingly availability of high-performance computers has made feasible detailed simulations of reactive transport in two- and three-dimensional domains (e.g. PFLOTRAN , retrieved on August 12, 2009), it is still not possible to re- solve practical problems with enough detail to capture all the scales of heterogeneity that are relevant for reactive transport (Steefel et al., 2005). Performance is also important when a large number of scenarios must be simulated. For example, because of our inability to observe all the scales of physical and chemical heterogeneity present in natural porous media, reactive transport simulations include a high degree of uncertainty. A standard way to deal with uncertainty is based on a Monte Carlo approach that involves simulating many equally probable scenarios to determine probability distributions for possible outcomes (Steefel et al., 2005). Simulating a large number of realizations is also a requirement of some methods to estimate the effects of variations in the input parameters on the model results (sensitivity analysis) or of some automatic parameter calibration frameworks (Hill and Tiedeman, 2007). Nowadays, high performance requires to use numerical methods that can be implemented in algorithms that are amenable to parallelization. While RWPT methods can be imple- mented using an embarrassingly parallel algorithm, efficient implementations of mesh- based algorithms, although possible, are more difficult to obtain (Mills et al., 2007). On the other hand, RWPT have very low convergence and very large numbers of particles (of the order of billions (Suciu et al., 2006)) are required to obtain accurate results, which counterbalances its parallel advantages. Finally, high-order multidimensional mesh-based solvers for advection are usually im- plemented using explicit schemes (Steefel and MacQuarrie, 1996). Explicit solvers have stability limits on the time step size given by the Courant–Friedrichs–Lewy (CFL) condi- tion, ∆t ≤ ∆/ |v|, where∆ is the cell size. Therefore, the stability limit is more restrictive for finer grids and finer discretizations require smaller time steps, hence, computational effort. Since mesh-based multidimensional solvers introduce a global coupling between concentration values, the stability limit is global and given by the maximum velocity in the grid. Thus, it is not possible to take advantage of the irregular velocity distribution to use greater time steps in slower areas of the domain. 16 1.4 Lagrangian Numerical Methods 1.4.1 Meshless Methods Kernel interpolation methods simulate mass transport using a collection of particles that move according to the velocity field and carry and exchange solute mass with surrounding particles. Particle locations are used as quadrature points to evaluate integral interpola- tions of variables and their derivatives. Importantly, these schemes are able to incorporate diffusive effects and mixing without using a grid or mesh, so they are also called mesh- less methods. Some examples of this type of method are Vortex methods (Cottet and Koumoutsakos, 2000) and particle strength exchange (PSE) method (Degond and Mas- Gallic, 1989a; Zimmermann et al., 2001). Methods based on this approach give accurate and stable results if a remeshing technique is used to control errors that result from ir- regular spatial particle distributions produced by non-uniform flow fields. The remeshing step introduces numerical dispersion that can be controlled but not completely avoided by using suitable interpolation schemes (Cottet and Koumoutsakos, 2000; Chaniotis et al., 2002). Smoothed Particle Hydrodynamics (SPH) methods are another type of kernel-based in- terpolation scheme (Gingold and Monaghan, 1977; Lucy, 1977). Cleary and Monaghan (1999) presented a SPH scheme that allows one to solve the multidimensional advection- dispersion equation, assuming isotropic dispersion, using an integral interpolation of the dispersion operator that it is supposed to be less sensitive to particle disorder than tra- ditional kernel interpolation schemes. Since the method can handle dispersion without remapping the concentration field onto a grid, it is free of numerical dispersion and grid orientation effects. As discussed in Chapter 2, a key property of kernel interpolation methods is that they track concentration values, in contrast to RWPT approaches which fundamentally track particles with fixed masses. This feature allows one to evaluate reactions at individual particles (e.g. Chaniotis et al., 2003). Heterogeneous reactions that include the solid (immobile) phase can be handled by introducing additional fixed particles (Tartakovsky et al., 2007). 17 1.4.2 Streamline-Based Simulations Streamline-based methods have been successfully used to simulate oil migration (Thiele et al., 1996, 1997) and multidimensional solute transport (Crane and Blunt, 1999; Di Do- nato and Blunt, 2004; Obi and Blunt, 2004, 2006). These methods use a numerical grid that adapts to the flow field, which reduces numerical dispersion and grid orientation effects. Because of its adaptation to the flow and its ability to minimize numerical dis- persion, the method is well suited for simulations of advection-dominated transport as found in heterogeneous porous media (Di Donato et al., 2003). Moreover, the use of streamlines allows the transformation of a multidimensional transport equation to a set of individual one-dimensional transport problems. The numerical solution of the result- ing set of one-dimensional transport problems allows the use of more efficient numerical solvers, more relaxed stability constraints and it is amenable to parallelization (Crane and Blunt, 1999; Bandilla et al., 2009). Because of the efficiency of the method, it is pos- sible to simulate large-scale domains with fine spatial and temporal resolution (Di Donato et al., 2003; Obi and Blunt, 2004, 2006). In addition, chemical reactions, including ho- mogeneous and heterogeneous reactions can be easily handled (Crane and Blunt, 1999; Di Donato and Blunt, 2004). 1.5 Discussion The previous analysis demonstrated that no single traditional numerical method presents the three main features sought in reactive transport simulations. Moreover, one of the conditions, monotonicity, is not satisfied by any of the current methods if anisotropic dispersion is considered. Because of the limitations of current numerical schemes, reactive transport simulations must make some trade-offs to obtain practical results. For example, low-order approx- imations for advection are preferred over more accurate high-order approximations, be- cause they do not suffer numerical oscillations. Similarly, simulations that include local- dispersion assume isotropic dispersion or remove the cross-terms to avoid introducing non-physical artifacts, e.g. negative concentrations. In both cases, monotonicity comes at the cost of tolerating additional numerical errors. In our opinion, most of the problems that affect current numerical methods are caused by the use of a grid or mesh, either to compute approximations of concentrations derivatives 18 as in finite difference or finite element methods or, to compute cell averaged concentration values as in RWPT, MOC and other hybrid approaches. As discussed, the use of a single multidimensional grid introduces grid orientations effects due to the non-uniform flow direction, artificial mixing because of the computation of cell averaged concentrations, and global stability restrictions. Therefore, it seems that to overcome many of the prob- lems that plague current numerical schemes, one should develop numerical methods that do not require a rigid multidimensional grid. Numerical methods based on a Lagrangian description of solute transport satisfy that condition and are attractive alternatives to simulate reactive solute transport. Because of their ability to control numerical dispersion and grid orientation effects and their efficiency, meshless methods and streamline-based simulations are attractive alter- natives for the simulation of reactive transport in natural porous media. However, they also present some deficiencies that can be problematic for their use in reactive transport codes. The standard SPH approximation for diffusion (Cleary and Monaghan, 1999) can only be used to simulate isotropic dispersion. Therefore, the use of this type of methods in reactive transport modeling in porous media requires deriving new expressions to simulate anisotropic dispersion. Although, approximations for longitudinal dispersion along individual streamlines are straightforward, transverse mixing between streamlines is more difficult to simulate. Two approaches have been used to incorporate transverse dispersion in streamline-based simu- lations. In the first one, solute transport is solved using a flow-oriented grid and transverse dispersion is included as a flux component perpendicular to the streamlines (Frind and Germain, 1986; Frind et al., 1987; Cirpka et al., 1999). This approach has been suc- cessfully used in two-dimensional simulations (Frind et al., 1987; Cirpka et al., 1999b), but it has not been extended to three-dimensions. A second alternative consists in using a hybrid approach (Obi and Blunt, 2004). First, advection is solved along streamlines. Then, concentration values are mapped onto a grid where a mesh-based solver is used to solve for dispersion. Finally, concentration values are interpolated back from the grid to the streamlines. The interpolation from and to streamlines introduces some numerical error that is difficult to quantify (Obi and Blunt, 2004). Because the interpolation must be done at each time step, the cumulative effect can be important even if an accurate interpolation scheme is used. The limitations of the meshless SPH and streamline-based methods for the simulation 19 of multidimensional reactive transport in porous media provide the motivation for the research in this thesis. 1.6 Objectives The main objective of this thesis is to develop, implement, and evaluate new numerical schemes based on meshless methods and streamline-based simulations to simulate reactive transport. The main objective includes the following specific objectives: 1. To derive and implement a meshless approximation for conservative transport in heterogeneous porous media. This includes deriving expressions to approximate isotropic and anisotropic local-scale dispersion. 2. To devise schemes to incorporate local-scale dispersion (longitudinal and transverse) in multidimensional streamline-based simulations. 3. To evaluate the newly derived numerical schemes in terms of accuracy, monotonicity and performance. 4. To compare the new schemes with others current numerical methods such as: high- order finite volume, method of characteristics and random-walk particle tracking methods. 5. To evaluate the suitability of using the new meshless and streamline-based schemes in reactive transport simulations. 1.7 Organization This thesis is organized in four additional chapters and four appendices. Chapters 2, 3, and 4 correspond to manuscripts that have been published or will be submitted for publication. Chapter 2 presents the application of a meshless numerical method based on smoothed particle hydrodynamics (SPH) for the simulation of conservative transport in heteroge- neous geological formations assuming isotropic dispersion. The chapter includes ana- lytical and numerical results that demonstrate that the new proposed scheme is stable, 20 accurate, and conserves global mass. Appendix A presents details of the derivation of the SPH-based numerical approximation. In Chapter 3, we extend the SPH-based approximation implemented in Chapter 2 to simulate anisotropic dispersion. In addition, we compare the new approximation with another meshless method (particle strength exchange) and a mesh-based finite volume scheme to simulate the dispersion of a two-dimensional contaminant plume under different scenarios. We conclude that, although attractive to simulate conservative transport, the new SPH-based approximation is unsuitable for reactive transport simulations because of spurious oscillations that arise if the dispersion tensor is anisotropic. The new numerical approximation is based on a SPH approximation for mixed second order derivatives, which is derived in detail in Appendix B. Chapter 4 presents the derivation of a new numerical scheme to incorporate dispersion – including transverse dispersion – in streamline simulations. A key element of the method is that dispersion is approximated in a flow oriented grid using a combination of a one-dimensional finite difference scheme and a meshless approximation for isotropic dispersion. We demontrate through analytical and numerical results that the resulting approximation is always monotonic and, hence, suitable for reactive transport simula- tions. Some key issues that arise in streamline-based simulations such as: streamline tracing, streamline spatial distribution and streamline discretization, are discussed in Appendix D. Finally, Chapter 5 summarizes the main conclusions of the three preceding chapters and includes recommendations for future research directions. 21 1.8 References Cleaning Up the Nation’s Wastes Sites: Markets and Technology Trends, Tech. Rep. 542R96005A, U.S. Environmental Protection Agency, 1996. Federal Contaminated Sites Inventory, http://www.tbs-sct.gc.ca/fcsi-rscf/home- accueil.aspx, retrieved on August 10, 2009. Ahlstrom, S., H. Foote, R. Arnett, C. Cole, and R. Serne, Multicomponent mass transport model: theory and numerical implementation (discrete-parcel-random-walk version), Tech. rep., BNWL-2127, Battelle Pacific Northwest Labs., Richland, Wash.(USA), 1977. Bandilla, K., A. Rabideau, and I. Janković, A parallel mesh-free contaminant transport model based on the Analytic Element and Streamline Methods, Adv. Water Resour., 32 , 1143–1153, 2009. Bear, J., Dynamics of fluids in porous media, Dover, 1988. Beckie, R., Scale dependence and scale invariance in hydrology, chap. Analysis of scale effects in large-scale solute-transport models, pp. 314–334, Cambridge University Press, 1998. Benson, D., and M. Meerschaert, Simulation of chemical reaction via particle tracking: Diffusion-limited versus thermodynamic rate-limited regimes, Water Resour. Res., 44 , 12, 2008. Boggs, J., S. Young, L. Beard, L. Gelhar, K. Rehfeldt, and E. Adams, Field study of dispersion in a heterogeneous aquifer 1. Overview and site description, Water Resour. Res., 28 , 3281–3291, 1992. Boris, J., and D. Book, Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works, J. Comput. Phys., 11 , 172–186, 1973. Burnett, R., and E. Frind, Simulation of contaminant transport in three dimensions: 2. Dimensionality effects, Water Resour. Res. WRERAQ, 23 , 695–705, 1987. Celia, M., T. Russell, I. Herrera, and R. Ewing, An Eulerian-Lagrangian localized adjoint method for the advection-diffusion equation, Water Resour., 13 , 187, 1990. 22 Chaniotis, A., D. Poulikakos, and P. Koumoutsakos, Remeshed smoothed particle hydro- dynamics for the simulation of viscous and heat conducting flows, J. Comput. Phys., 182 , 67–90, 2002. Chaniotis, A. K., C. E. Frouzakis, J. C. Lee, A. G. Tomboulides, and K. Poulikakos, D. AU Boulouchos, Remeshed smoothed particle hydrodynamics for the simulation of laminar chemically reactive flows, J. Comput. Phys., 191 , 1–17, 2003. Cirpka, O., E. Frind, and R. Helmig, Numerical methods for reactive transport on rect- angular and streamline-oriented grids., Adv. Water Res., 22 , 711–728, 1999a. Cirpka, O., E. Frind, and R. Helmig, Numerical simulation of biodegradation controlled by transverse mixing, J. Contam. Hydrol., 40 , 159–182, 1999b. Cirpka, O. A., A. Olsson, Q. S. Ju, M. A. Rahman, and P. Grathwohl, Determination of transverse dispersion coefficients from reactive plume lengths, Ground Water , 44 , 212–221, 2006. Cleary, P. W., and J. J. Monaghan, Conduction modelling using smoothed particle hy- drodynamics, J. Comput. Phys., 148 , 227–264, 1999. Cottet, G., and P. Koumoutsakos, Vortex methods: Theory and practice., Cambridge University Press, 2000. Cox, R., and T. Nishikawa, A new Total Variation Dimishing scheme for the solution of advective-dominant solute transport, Water Resour. Res., 27 , 2645–2654, 1991. Crane, M., and M. Blunt, Streamline-based simulation of solute transport,Water Resour. Res., 35 , 3061–3078, 1999. De Simoni, M., J. Carrera, X. Sánchez-Vila, and A. Guadagnini, A procedure for the solu- tion of multicomponent reactive transport problems,Water Resour. Res., 41 , W11,410, 2005. Degond, P., and S. Mas-Gallic, The weighted particle method for convection-diffusion equations. Part 1: The case of an isotropic viscosity, Math. Comput., 53 , 485–507, 1989a. Degond, P., and S. Mas-Gallic, The weighted particle method for convection-diffusion equations. II: The anisotropic case, Math. Comp, 53 , 485,508, 1989b. 23 Delay, F., P. Ackerer, and C. Danquigny, Simulating Solute Transport in Porous or Fractured Formations Using Random Walk Particle Tracking: A Review, Vadose Zone J., 4 , 360–379, 2005. Di Donato, G., and M. Blunt, Streamline-based dual-porosity simulation of reactive trans- port and flow in fractured reservoirs, Water Resour. Res., 40 , 2004. Di Donato, G., E. Obi, and M. Blunt, Anomalous transport in heterogeneous media demonstrated by streamline-based simulation, Geophys. Res. Lett., 30 , 1608, 2003. Edwards, M., and H. Zheng, A quasi-positive family of continuous Darcy-flux finite- volume schemes with full pressure support, J. Comput. Phys., 2008. Freyberg, D., A natural gradient experiment on solute transport in a sand aquifer: 2. Spa- tial moments and the advection and dispersion of nonreactive tracers, Water Resour. Res, 22 , 2031–2046, 1986. Frind, E., and D. Germain, Simulation of contaminant plumes with large dispersive con- trast: Evaluation of alternating direction galerkin models, Water Resour. Res. WR- ERAQ, 22 , 1986. Frind, E., E. Sudicky, and S. Schellenberg, Micro-scale modelling in the study of plume evolution in heterogeneous media, Stoch. Hydrol. Hydraul., 1 , 263–279, 1987. Gingold, R. A., and J. J. Monaghan, Smoothed particle hydrodynamics: Theory and application to non-spherical stars, Mon. Not. R. Astron. Soc., 181 , 375–389, 1977. Ham, P., R. Schottinga, H. Prommerb, and G. Davisc, Effects of hydrodynamic dispersion on plume lengths for instantaneous bimolecular reactions, Adv. Water Resour., 27 , 803–813, 2004. Harten, A., and P. Lax, On a class of high resolution total-variation-stable finite-difference schemes, SIAM J. Numer. Anal., pp. 1–23, 1984. Herrera, P., and A. Valocchi, Positive solution of two-dimensional solute transport in heterogeneous aquifers, Ground Water , 44 , 803–813, 2006. Hess, K., J. Davis, D. Kent, and J. Coston, Multispecies reactive tracer test in an aquifer with spatially variable chemical conditions, Cape Cod, Massachusetts: Dis- persive transport of bromide and nickel, Water Resour. Res., 38 , 1161–1177, 2002. 24 Hill, M., and C. Tiedeman, Effective groundwater model calibration: With analysis of data, sensitivities, predictions, and uncertainty, Wiley-Interscience, 2007. Jose, S. C., M. A. Rahman, and O. A. Cirpka, Large-scale sandbox experiment on lon- gitudinal effective dispersion in heterogeneous porous media, Water Resour. Res., 40 , W12,415, 2004. Kitanidis, P. K., The concept of the dilution index, Water Resour. Res., 30 , 2011–2026, 1994. Konikow, L., D. Goode, G. Hornberger, and G. Survey, A Three-dimensional Method-of- characteristics Solute-transport Model (MOC3D), US Geological Survey, 1996. Le Potier, C., Finite volume monotone scheme for highly anisotropic diffusion operators on unstructured triangular meshes, Comptes Rendus Mathématique, 341 , 787–792, 2005. Lichtner, P., S. Kelkar, and B. Robinson, New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking., Water Resour. Res., 38 , 1146, 2002. Lipnikov, K., D. Svyatskiy, and Y. Vassilevski, Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes, J. Comput. Phys., 228 , 703–716, 2009. Lucy, L., A numerical approach to the testing of the fission hypothesis, Astron. J., 82 , 1013–1024, 1977. Mackay, D., D. Freyberg, P. Roberts, and J. Cherry, Natural Gradient Experiment on Solute Transport in a Sand Aquifer: 1. Approach and Overview of Plume Movement, Water Resour. Res. WRERAQ, 22 , 1986. Mayer, K., E. Frind, and D. Blowes, Multicomponent reactive transport modeling in vari- ably saturated porous media using a generalized formulation for kinetically controlled reactions, Water Resour. Res., 38 , 1174, 2002. Mills, R., C. Lu, P. Lichtner, and G. Hammond, Simulating subsurface flow and transport on ultrascale computers using PFLOTRAN, in Journal of Physics: Conference Series, vol. 78, p. 012051, Institute of Physics Publishing, 2007. 25 Mlacnik, M., and L. Durlofsky, Unstructured grid optimization for improved monotonicity of discrete solutions of elliptic equations with highly anisotropic coefficients, J. Comput. Phys., 216 , 337–361, 2006. Molz, F., and M. Widdowson, Internal inconsistencies in dispersion-dominated models that incorporate chemical and microbial kinetics, Water Resour. Res., 24 , 1988. Neuman, S., A Eulerian-Lagrangian numerical scheme for the dispersion-convection equa- tion using conjugate space-time grids, J. Comput. Phys., 41 , 1981. Neuman, S., Adaptive Eulerian-Lagrangian finite element method for advection- dispersion, Int. J. Numer. Meth. Engng., 20 , 321–37, 1984. Nordbotten, J., and I. Aavatsmark, Monotonicity conditions for control volume methods on uniform parallelogram grids in homogeneous media, Computat. Geosci., 9 , 61–72, 2005. Obi, E., and M. Blunt, Streamline-based simulation of advective-dispersive solute trans- port, Adv. Water Resour., 27 , 913–924, 2004. Obi, E. I., and M. J. Blunt, Streamline-based simulation of carbon dioxide storage in a North Sea aquifer, Water Resour. Res., 42 , W03,414, 2006. Oya, S., and A. J. Valocchi, Transport and biodegradation of solutes in stratified aquifers under enhanced in situ bioremediation conditions, Water Resour. Res., 34 , 3323–3334, 1998. PFLOTRAN, Scaling - application parallel performance, http://ees.lanl.gov/source/orgs/ees/pflotran/simscaling.shtml, retrieved on August 12, 2009. Pickens, J., and G. Grisak, Scale-dependent dispersion in stratified granular aquifer., Water Resour. Res., 17 , 1191–1211, 1981. Pruess, K., TOUGH2: A general-purpose numerical simulator for multiphase fluid and heat flow, Tech. Rep. LBL-29400 , Lawrence Berkeley National Laboratory, 1991. Rubin, J., Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions, Water Resour. Res., 19 , 1983. 26 Russell, T., and M. Celia, An overview of research on Eulerian-Lagrangian localized adjoint methods (ELLAM), Adv. Water Resour., 25 , 1215–1231, 2002. Salamon, P., D. Fernàndez-Garcia, and J. Gómez-Hernández, A review and numerical assessment of the random walk particle tracking method., J. Contam. Hydrol., 87 , 277–305, 2006. Steefel, C., and K. MacQuarrie, Approaches to modeling of reactive transport in porous media, Reviews in Mineralogy and Geochemistry, 34 , 85–129, 1996. Steefel, C., and K. Maher, Fluid-rock interaction: A reactive transport approach, Tech. Rep. LBNL1798E , Lawrence Berkeley National Laboratory, 2009. Steefel, C., D. DePaolo, and P. Lichtner, Reactive transport modeling: An essential tool and a new research approach for the Earth sciences, Earth Planet. Sci. Lett., 240 , 539–558, 2005. Suciu, N., C. Vamos, J. Vanderborght, H. Hardelauf, and H. Vereecken, Numerical in- vestigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res, 42 , 1–17, 2006. Sudicky, E., A natural gradient experiment on solute transport in a sand aquifer: Spa- tial variability of hydraulic conductivity and its role in the dispersion process, Water Resour. Res., 22 , 2069–2082, 1986. Tartakovsky, A., P. Meakin, T. Scheibe, and B. Wood, A smoothed particle hydrody- namics model for reactive transport and mineral precipitation in porous and fractured porous media, Water Resour. Res., 43 , W05,437, 2007. Thiele, M., R. Batycky, and M. Blunt, Simulating flow in heteroneous systems using streamtube and streamlines, SPE Reservoir Engineering, pp. 5–12, 1996. Thiele, M., R. Batycky, and M. Blunt, A streamline-based 3D field-scale compositional reservoir simulator, Soc. Petrol. Eng. J., 1997. Tompson, A., Numerical simulation of chemical migration in physically and chemically heterogeneous porous media, Water Resour. Res., 29 , 3709–3726, 1993. Tompson, A., and L. Gelhar, Numerical simulation of solute transport in randomly het- erogeneous porous media, Water Resour. Res., 26 , 2541–2562, 1990. 27 Valocchi, A., Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils, Water Resour. Res., 21 , 1985. Weeks, S., and G. Sposito, Mixing and stretching efficiency in steady and unsteady groundwater flows, Water Resour. Res., 34 , 3315–3322, 1998. White, M., M. Oostrom, and R. Lenhard, Modeling fluid flow and transport in variably saturated porous media with the STOMP simulator. 1. Nonvolatile three-phase model description, Adv. Water Resour., 18 , 353–364, 1995. Yee, H., R. Warming, and A. Harten, Implicit total variation diminishing (TVD) schemes for steady-state calculations, J. Comput. Phys, 57 , 327–360, 1985. Yeh, G., and V. Tripathi, A critical evaluation of recent developments in hydrogeochem- ical transport models of reactive multichemical components, Water Resour. Res., 25 , 93–108, 1989. Yuan, G., and Z. Sheng, Monotone finite volume schemes for diffusion equations on polygonal meshes, J. Comput. Phys., 2008. Zalesak, S., Fully multidimensional flux-corrected transport algorithms for fluids, J. Com- put. Phys., 31 , 335–362, 1979. Zheng, C., and S. M. Gorelick, Analysis of solute transport in flow fields influenced by preferential flowpaths at the decimeter scale, Ground Water , 41 , 142–155, 2003. Zheng, C., and P. Wang, MT3DMS: A Modular Three-Dimensional Multispecies Trans- port Model for Simulation of Advection, Dispersion, and Chemical Reactions of Con- taminants in Groundwater Systems; Documentation and User’s Guide, Contract Report SERDP-99-1, US Army Engineer Research and Development Center, Vicksburg, MS , 1999. Zimmermann, S., P. Koumoutsakos, and W. Kinzelbach, Simulation of pollutant trans- port using a particle method, J. Comput. Phys., 173 , 322–347, 2001. 28 Chapter 2 A Meshless Method to Simulate Solute Transport in Heterogeneous Porous Media1 2.1 Introduction 2.1.1 Background Contaminant transport in natural aquifers is a complex, multiscale process that is fre- quently studied using numerical methods. Conservative solute transport is typically modeled using the advection-dispersion equation (ADE). Despite the large number of available numerical methods that have been developed to solve it, the accurate numerical solution of the ADE still presents formidable challenges. In particular, current numeri- cal solutions of multidimensional advection-dominated transport in non-uniform velocity fields are affected by one or all of the following problems: numerical dispersion that in- troduces artificial mixing and dilution, grid orientation effects, and unphysical numerical oscillations (Herrera and Valocchi, 2006). To correctly capture the basic mechanisms that control conservative solute transport in natural aquifers an ideal numerical method should be able to: (i) accurately capture 1A version of this chapter has been published. P. Herrera, M. Massabo, and R. Beckie (2009) A meshless method to simulate solute transport in heterogeneous porous media. Adv. Water Resour., 32:413–429. 29 the effect of small-scale velocity fluctuations upon solute distribution, (ii) simulate the effect of small values of local-scale dispersion in advection-dominated transport, (iii) reproduce small-scale concentration fluctuations and (iv) allow the efficient simulation of problems at moderate to large scales. In what follows we only discuss numerical methods to simulate conservative solute transport but the same set of requirements should be satisfied by any successful numerical method used to simulate reactive transport. The objective of this paper is to develop and test a meshless method to simulate contam- inant transport in porous media. This work was primarily motivated by our experience in theoretical investigations of solute mixing and plume dilution in heterogeneous porous media. Those investigations require an efficient numerical method that is able to accu- rately simulate solute transport in multidimensional non-uniform velocity fields. We are also interested in developing a meshless method that can be used to simulate conservative and non-conservative solute transport with minimum modification. The main contributions of this paper are: (i) to derive a meshless approximation for the dispersion operator in heterogeneous porous media flow, (ii) to compare the meshless approximation with other numerical methods traditionally used to simulate conserva- tive solute transport in heterogeneous multidimensional velocity fields, (iii) to study the convergence properties of the meshless approximation of the dispersion operator for different spatial node distributions, and (iv) to demonstrate that the proposed mesh- less approximation can be used to solve the ADE in different scenarios ranging from advection-dominated to dispersion-dominated solute transport with accuracy better than or comparable to standard numerical methods. Although we address only conservative transport in this paper we highlight the advantages of the proposed scheme for simulating reactive transport. 2.1.2 Numerical Methods A detailed discussion of every numerical method used to solved the ADE is beyond the scope of this manuscript, however we briefly discuss some of them to motivate the development of our new meshless method. Grid or mesh-based methods such as finite difference, finite volume or finite element meth- ods are relatively easy to implement, their convergence, stability and accuracy properties are well understood, and it is possible to develop formulations that are mass conserva- tive. On the other hand, mesh-based methods have difficulty in accurately simulating 30 multidimensional advection-dominated problems because of numerical dispersion that in- troduces excessive artificial mixing and dilution of the contaminant plume. Therefore, grid-based methods are only advised for problems with low Péclet number (Zheng and Wang, 1999). Hybrid schemes were developed to address the limitations of grid-based methods. The following hybrid schemes are based on the same general concept – the use of streamlines or particles to handle advection and a grid-based method to handle dispersion: hybrid Eulerian-Lagrangian methods (Neuman, 1981, 1984), method of characteristics (MOC) and hybrid method of characteristics (HMOC) in the popular MT3DMS (Zheng and Wang, 1999) and MOC3D (Konikow et al., 1996), Eulerian-Lagrangian localized adjoint methods (ELLAM) (Celia et al., 1990; Russell and Celia, 2002), and three-dimensional hybrid streamline-grid approaches (Obi and Blunt, 2004). At each time step, these meth- ods solve the advection–dispersion equation in two steps. First, changes in concentrations due to advection are computed using some suitable scheme such as particle-tracking or by solving the transport equation along streamlines. Then, concentration values are in- terpolated onto a grid. Next, grid concentration values are used to solve for dispersion, and eventually reactions, using some traditional mesh-based solver. Multiple variations of this approach exist depending on the interpolation methods and tracking algorithm. However, all of them require mapping concentrations between cell centers and particle locations or streamline nodes. This remapping step introduces numerical dispersion that is difficult to quantify and control in multidimensional simulations. Moreover, the nu- merical dispersion due to the remapping step is more important in simulations with large grid Péclet number where the effect of dispersion is not enough to smooth the sub-grid scale concentration distribution. The limitations of grid-based and hybrid methods to simulate advection dominated prob- lems have motivated the development of Lagrangian and meshless methods including those based on a random-walk and those based on integral interpolations. Random-walk particle-tracking methods (RWPT) have long been used to simulate con- servative solute transport in porous media (Ahlstrom et al., 1977; Smith and Schwartz , 1980; Pickens and Grisak, 1981). In this type of model, solute mass is distributed among a set of particles that move carried by the flow velocity and by a random drift that mod- els dispersive transport (Delay et al., 2005; Salamon et al., 2006). Solute concentrations are estimated by averaging the mass contained in the particles found in some specified volume. Therefore, concentration values depend upon the total number of particles, size of the averaging volume, and spatial particle distribution. The numerical precision of 31 the computed concentration is limited by the finite number of particles used in a simula- tion and the calculated concentration usually exhibits numerical oscillations that can be amplified in presence of non-linear reactions (Tompson and Dougherty, 1992; Tompson, 1993). Therefore, the use of the method is limited to conservative transport or to reac- tive transport simulations with simple reactions that can be modeled by changing the state or phase of individual particles, e.g. sorption (Valocchi and Quinodoz , 1989). The popularity of RWPT is due to its natural capacity to accurately simulate advection, ease of implementation, and relatively moderate computational requirements. Because of its advantages RWPT has become the de facto standard method used in numerical studies of plume spreading and dilution. In contrast to RWPT methods that simulate dispersion in a fluid by a random movement of particles that carry solute mass, kernel interpolation methods simulate mass transport using a collection of particles that move according to the velocity field and carry and ex- change solute mass with surrounding particles. Particle locations are used as quadrature points to evaluate integral interpolations of variables and their derivatives. Importantly, these schemes are able to incorporate diffusive effects and mixing without using a grid or mesh, so they are also called meshless methods. Some examples of this type of method are vortex methods (Cottet and Koumoutsakos, 2000) and particle strength exchange (PSE) method (Degond and Mas-Gallic, 1989a). Zimmermann et al. (2001) present, to the best of our knowledge, the only application of this type of method to simulate solute transport in porous media at the continuum scale. Their results indicate that the method gives accurate and stable results for the flow configuration considered if a remeshing tech- nique is used to control errors that result from irregular spatial particle distributions. As in other methods, the remeshing step introduces numerical dispersion that can be con- trolled but not completely avoided by using a suitable interpolation scheme (Cottet and Koumoutsakos, 2000; Chaniotis et al., 2002). Smoothed particle hydrodynamics (SPH) methods are another type of kernel-based in- terpolation scheme (Gingold and Monaghan, 1977; Lucy, 1977). Cleary and Monaghan (1999) presented a SPH scheme that allows one to solve the multidimensional ADE us- ing an integral interpolation of the dispersion operator that is less sensitive to particle disorder than traditional kernel interpolation schemes. Since the method can handle dispersion without remapping the concentration field onto a grid, it is free of numerical dispersion and grid orientation effects. As we show later, a key property of SPH is that the method tracks concentration values, in contrast to RWPT approaches which funda- mentally track particles with fixed masses. This feature of SPH allows one to simulate 32 reactive transport directly with individual particles (e.g. Chaniotis et al., 2003). We next propose a method to simulate solute transport in porous media based on the SPH formalism. We show that this SPH-based method is advantageous because of its inherent ability to solve advection, its capacity to approximate dispersion in a meshless fashion, and its ability to reproduce smooth fine-scale concentration distributions appropriate for reactive transport simulations. After deriving the method in Section 2.2, we compare and contrast it to a conventional RWPT method and other traditional mesh-based methods in Section 2.3. 2.2 Monte Carlo SPH method In application of meshless methods such as SPH it is appropriate to reformulate the ADE in terms of a Lagrangian coordinate system as dr dt = v (r) (2.1) dC dt r = ∇ · (D∇C)|r (2.2) where C is the solute concentration, D is the local-scale dispersion tensor, v is the pore water velocity, and r is the position vector of a small fluid volume or material point. As we describe next, in SPH methods the concentration field is represented using a set of particles that carry concentration information and are distributed through the domain – even in areas where solute concentration is zero. Practically all meshless methods including SPH use a particle-tracking approach to integrate (2.1) in the same way as done in RWPT (Delay et al., 2005; Salamon et al., 2006). The key distinction between SPH and RWPT is how dispersion is approximated. In SPH methods, dispersion – the solution of equation (2.2) – is evaluated using a kernel interpolation approximation. The “smoothed” in SPH comes from the representation of a scalar or vector field by a smoothed integral interpolation. The smoothed interpolation AS(r) of a field A(r) is defined as the integral (Gingold and Monaghan, 1977) AS(r) = ˆ A(r)W (r− r, h) dr (2.3) 33 where W (r− r, h) is a kernel function with compact support around r and smoothing length h that satisfies some normalization properties (Gingold and Monaghan, 1977). Spline polynomials with finite support are usually used as kernel functions because of their practical advantages (Price, 2004). For those functions the compact support volume of the kernel depends upon h. Figure 2.1 shows values of the Gaussian and cubic-spline kernels and their derivatives as function of the smoothing length. −3 −2 −1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 r/h W (r/ h) −3 −2 −1 0 1 2 3 −0.8 −0.6 −0.4 −0.2 0.0 r/h W ’(r /h) −3 −2 −1 0 1 2 3 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 r/h W ’’( r/h ) −3 −2 −1 0 1 2 3 −2.0 −1.5 −1.0 −0.5 0.0 r/h F( r/h ) Figure 2.1: Kernel function W, first derivative W’, second derivative W’’ and sym- metric F function defined in Appendix A, for Gaussian (solid line) and cubic spline (dashed line) kernels. The numerical approximation of AS can be evaluated using a Monte Carlo integration scheme by sampling the integrand A(r)W (r) at a limited set of disordered points or particle locations. To evaluate the integral one must consider that the set of scattered points is not uniformly distributed, therefore their spatial distribution must be explicitly taken into account in the integral evaluation to get the following modified form of (2.3), As(r) = ˆ A(r) p(r)W (r− r ) p(r)dr (2.4) 34 where p(r) is the probability density of finding a particle in a given unit volume with units of one over volume [1/L3] (Gingold and Monaghan, 1982; Tartakovsky and Meakin, 2005). Thus, p(r)dr can be interpreted as a non-uniform density to sample the modified integrand A(r)W (r)/p(r) (Press et al., 1992, p. 316). The exact evaluation of p(r) for a set of scattered points in multiple dimensions is a very difficult task, but it can be estimated using a density estimation by kernels approach (Schaback and Wendland, 2006), which yields, p(r) = 1 Np Np j=1 W (r− rj) (2.5) Finally, the Monte Carlo approximation of (2.4) is AN(r) = 1 Np Np j A(rj)W (r− rj) p(rj) ±O 1 Np (2.6) where Np is the total number of points that effectively contribute to the integral. This last expression is a valid approximation of the integral for any set of scattered points. We note that in practice the error estimate in (2.6) is an upper bound and that the actual error also depends upon p(r), i.e. the spatial distribution of the points, and the local smoothness properties of A. For example, numerical studies have shown that the actual error is much smaller for points that are reasonably well distributed (Cleary and Monaghan, 1999; Monaghan, 2005). Traditional SPH simulations used to solve hydrodynamics equations consider that the set of locations rj represent the positions of a set of fluid particles with constant mass mj. In that case, the fluid mass density ρ(r) is proportional to the particle density p(r) and it can be estimated by using an expression similar to (2.5) where the mass of individual particles appears explicitly (Monaghan, 1992, 2005). We note that the formulation given by (2.5) and (2.6) is equivalent to the standard SPH formulation if the total mass in the system is one and it is uniformly distributed among particles. This is consistent with SPH simulations where, in general, the fluid mass assigned to each particle is considered as a constant scaling parameter that is set at the beginning of individual simulations to tune the spatial fluid density distribution (Dilts, 1999). 35 There are two SPH approaches to approximate transport equations involving second order derivatives. First, second order derivatives can be easily evaluated by differentiating (2.6). The resulting expression involves the second derivative of the kernel, so it is very sensitive to particle disorder (Cleary and Monaghan, 1999). In that case particle positions must be periodically reinitialized on a regular lattice to achieve acceptable accuracy (Chaniotis et al., 2002). A second approach is based on an integral approximation of the dispersive fluxes that depends only upon the first derivative of the kernel (Brookshaw, 1985; Cleary and Monaghan, 1999). Considering isotropic dispersion, i.e. D(r) = D(r)I, where I is the identity matrix, the approximation of (2.2) is dCi dt = 1 Np Np j=1 1 pj (Dj +Di) (Cj − Ci)F (rj − ri) (2.7) where Ci is the solute concentration at position ri and F (rj − ri) is a function of the separation vector and first derivative of the kernel that has spherical symmetry. Appendix A presents the derivation of (2.7). Equation (2.7) indicates that the magnitude of the contribution of solute from particle j to particle i is equal to the contribution of particle j to particle i, i.e. the mass flux is anti- symmetric fij = −fji, only if pi = pj . In general, for a set of irregularly spaced particles pi = pj , so (2.7) does not satisfy a basic property of the dispersion operator (Kuzmin and Turek, 2002). Thus, it is necessary to replace the denominator by a symmetric approximation of the form p̂ij = g(pi, pj), e.g. arithmetic or harmonic average. In general, this correction is relatively small because particles that effectively contribute to the summation are within few smoothing lengths due to the compact support property of the kernel. Similar corrections to recover symmetry are used in standard SPH simulations that consider variable smoothing lengths (Monaghan, 2005). Our MC-SPH method is based upon this formulation to approximate dispersive fluxes that results in dCi dt = 1 Np Np j=1 1 p̂ij (Dj +Di) (Cj − Ci)F (rj − ri) (2.8) We refer to this approximation as the Monte Carlo SPH (MC-SPH) formulation for dis- persion. Equation (2.8) satisfies two important physical constraints. First, dispersive fluxes are anti-symmetric. Second, mass transfer occurs from higher to lower concentra- tions because for typical kernels, F (rj − ri) < 0 as shown in Figure 2.1. 36 In standard SPH simulations there is another approach to recover the symmetry of the fluxes that consists in including the fluid density on the right hand side of (2.2). In this case, density is placed inside operators following the “second golden rule of SPH” (Monaghan, 1992). The resulting numerical approximation is dCi dt = Np j=1 mj ρiρj (ρjDj + ρiDi) (Cj − Ci)F (rj − ri) (2.9) where ρi is the numerical approximation of the fluid mass density at position ri. Equation (2.9) is the standard SPH approximation for dissipative or dispersive transport and it has been used to simulate heat conduction in compressible gases (Cleary and Monaghan, 1999; Español and Revenga, 2003; Jubelgas et al., 2004); to simulate viscous effects in low Reynolds number flows (Morris et al., 1997; Zhu et al., 1999); and to simulate solute dispersion at the pore scale (Zhu and Fox , 2002; Tartakovsky and Meakin, 2005). In what follows we refer to this approximation as weakly compressible SPH (WC-SPH) formulation for dispersion. Although this approach is reasonable in simulations that consider variable fluid density so that density must be explicitly incorporated inside the dispersion operator, its use in incompressible flow simulations with constant density is, at least, questionable. Moreover, as shown by equation (2.8), it is not necessary to use such a pragmatic approach to develop a numerical formulation that satisfies basic physical principles and conserves solute mass. We close this section by summarizing the main distinctive properties of the MC-SPH formulation. The approach considers that fluid particles represent a constant fluid volume larger than a representative elementary volume (REV) such that Darcy’s velocities can be defined and that particle trajectories can be computed by integrating (2.1). Solute mass is distributed among a set of particles that carry concentration values. Local-scale dispersion that occurs at scales much smaller than the REV, is modeled as a Fickian solute mass transfer process between particles. Numerically, local-dispersion is approximated by a local integral interpolation of the dispersion operator in (2.2). We note that the particle fluid volume does not explicitly appear in the numerical formulation and that, from a practical point of view, fluid particles can be regarded as numerical nodes. 2.2.1 Time Integration The time integration of the system of equations (2.1) and (2.2) requires the use of a sequential procedure. First, at the beginning of each time step node locations and con- 37 centrations are recorded. Then, new locations are calculated by integrating (2.1) using an explicit time-marching scheme (e.g. Runge–Kutta methods) or a particle-tracking al- gorithm (e.g. Pollock, 1988). After the new locations are computed, new concentrations are calculated by integrating (2.2). This term can be integrated using explicit or implicit schemes. For example, using a first-order approximation it can be approximated by Ct+∆ti − Cti ∆t = j αij C∗j − C∗i (2.10) where αij = 1 p̂j (Di +Dj)F (rj − ri) ≥ 0 (2.11) where C∗i = Cti or C∗i = Ct+∆ti for explicit and implicit time integration, respectively. Then, the first-order implicit approximation is given by 1 +∆t j αij Ct+∆ti − ∆t j αij Ct+∆tj = Cti (2.12) It is easy to demonstrate that this integration scheme is unconditionally stable and posi- tivity preserving. Although possible, implicit schemes are seldom used because of compu- tational overhead. Since nodes move with the flow, the connectivity list, i.e. the number of nodes within the kernel support volume of a given node changes at each time step. The memory requirements to store the associated matrix and the computational cost to generate it and computing its inverse can be quite large, depending upon the average number of nodes per smoothing length. That is why we used conditionally stable ex- plicit time integration schemes in the simulations presented in Section 2.3. We motivate the discussion about the stability of explicit schemes by writing the first-order explicit approximation of (2.10) Cti = 1−∆t j αij Cti + ∆t j αij Ctj (2.13) which is stable and positivity preserving if 1/αij ≥ ∆t. Empirical tests have shown that other explicit solutions are stable for time steps∆t such that (Cleary and Monaghan, 1999) 38 ∆t ≤ h 2 D (2.14) where is a coefficient that depends upon the kernel function. In our simulations we used a cubic-spline kernel and = 0.1 to get stable results using an explicit second-order Runge–Kutta scheme. From a physical point of view, equation (2.12) indicates that the time step must be smaller than a dispersion time scale given by the dispersion coefficient and the kernel support volume. 2.2.2 Accuracy and Spatial Resolution Errors in meshless approximations based on kernel interpolants come from two sources (Brackbill, 2005; Quinlan et al., 2006). The integral interpolation in (2.3) introduces a smoothing error that depends on the shape and smoothing length of the kernel. Addition- ally, there is a numerical integration error that depends upon the number and location of the nodes and the smoothness of the real function (Schaback and Wendland, 2006). There is a tradeoff between both sources of error because the accuracy of the smoothed quantity increases as the smoothing length decreases, while the numerical integration er- ror increases as the number of nodes per support volume decreases. Therefore, the only way to simultaneously reduce both sources of error is to decrease the smoothing length while increasing the total number of nodes to keep the same average number of nodes per support volume. In general, for a given function and node distribution there is a critical smoothing length hc such that for h > hc the smoothing error dominates and that for h < hc the integration error is more important. The determination of hc for irregularly spaced points is very difficult, if not impossible. The numerical integration error also depends upon the regularity of the node distribution. Meshless approximations such as SPH methods produce very accurate results in situa- tions where particles are regular or uniformly distributed (Monaghan, 2005). In those situations the error affecting the simulation is much smaller than the theoretical error es- timate of (2.6) which considers a random particle distribution (Monaghan, 2005). When particles are uniformly distributed the leading error term is due to the interpolation er- ror and is controlled by the kernel smoothing length that sets the spatial resolution. For example, numerical studies have shown that for a set of equispaced particles the error for (2.9) converges with h2 (Cleary and Monaghan, 1999). 39 The previous discussion indicates that the accuracy of the MC-SPH method to simulate local-dispersion will evolve during the simulation. At early times the node distribution is similar to the initial regular distribution, so the leading error term is due to the interpolation approximation. As nodes are redistributed in space by the non-uniform flow velocity, they become clustered in different zones, so the numerical integration error becomes much more important. Then, the accuracy of the method during the simulation must be controlled by an appropriate choice of the kernel smoothing length and the number of nodes per kernel support volume or initial average spacing. The choice of those parameters is not trivial, particularly in multidimensional simulations, and requires some trial and error. In the simulations presented below, we selected those parameters using the following steps: (i) setup an initial node configuration, (ii) simulate advection only, (iii) check node distribution and, particularly, number of nodes per kernel support volume, (iv) if node distribution was considered too sparse, increase the initial number of nodes, (v) repeat until an acceptable final node distribution is produced. Because particle tracking is very efficient, the determination of the optimum initial node configuration demanded little time and effort relative to the overall simulation time. 2.2.3 Mass Conservation In the proposed MC-SPH scheme solute mass is distributed in space as a finite set of concentration values at node positions, thus the total mass in the system equal to the integral of the concentration over the domain can be approximated as M = ˆ Ω C(r)dr ≈ CV +O 1 Np (2.15) where we have used a Monte Carlo integration approach and V is the volume of the domain Ω. As discussed above, the accuracy of the integral depends upon the total number of points Np and their spatial distribution. In most practical SPH simulations the number of particles is quite large and (2.15) is a good approximation. Therefore, we can study the evolution of the total solute mass in the system by writing dM dt ≈ V d C dt = V Np Np i=1 dCi dt (2.16) 40 Finally, substituting expression (2.10) used to compute the temporal derivative of the concentration, we get dM dt = V Np i j αij (Cj − Ci) = 0 (2.17) since αij = −αji, and thus solute mass is globally conserved. In the analyses that follow we use this expression to characterize mass balance. 2.3 Numerical evaluation of the MC-SPH method To test the capacity of the MC-SPH method to provide reasonably accurate solutions for dispersive transport we compare it with the analytical solutions of simple one- and two-dimensional dispersion problems and with other numerical solutions for simulating advective-dispersive transport in non-uniform velocity fields. 2.3.1 One-Dimensional Dispersion We consider a simple one-dimensional problem to illustrate the behavior of the error affecting our new MC-SPH approximation in (2.8) as function of the particle distribution and smoothing length. We simulate the dispersion of a one-dimensional Gaussian plume where the concentration as function of position and time is given by C(x, t) = C0s0 s e −(x−x0)2 2s2 (2.18) where s = s20 + 2Dt, C0 is the maximum initial concentration, s0 is a constant that controls the size of the initial plume, and x0 is the position of the plume center. Particles were initially distributed over a regular equispaced grid with spacing ∆x. To study the effect of the particle distribution, we generated a non-uniform particle distri- bution by adding a normally distributed perturbation with standard deviation σ. Then, we computed the numerical solution using a Gaussian kernel with cutoff at 4h, i.e. only particles within 4h contribute to the kernel summation. The error introduced by this approximation is small because of the rapid falloff of the Gaussian function. 41 Figure (2.2) shows the maximum normalized error defined as error = CNumerical − CAnalytical2 /N (2.19) where N = Np is the number of particles, versus h for two different ratios ∆x/h for solutions computed using our MC-SPH and the traditional WC-SPH, formulations. For uniform particle distribution (i.e., σ/∆x = 0) and ∆x/h = 0.66 the error increases as h2. As particles become disordered the error increases for any value of h. Particle disorder is less important for large smoothing lengths and errors for different σ are similar. Larger values of the ratio ∆x/h, which is equivalent to fewer particles per support volume, produces larger error even for the uniform particle distribution. For large particle disorder (σ/∆x = 1) the interpolation error is dominant and the error of the numerical solutions is almost independent of h. Figure (2.2) also shows that the use of the new MC-SPH approximation instead of the standard WC-SPH does not make a difference as both solutions produce results with similar accuracy. 2.3.2 Two-Dimensional Dispersion In this section we consider the simulation of the dispersion of a two-dimensional Gaus- sian plume using the RWPT and MC-SPH numerical methods. Despite the fact that both methods are based on a Lagrangian formulation of the solute transport problem and use a particle-tracking algorithm to integrate the advection equation, there are im- portant differences in their conceptual approaches, accuracy, numerical implementation, and computational performance that must be considered to evaluate their merits. Table 2.1 summarizes the main differences between both methods and Appendix C gives details about our implementation of the RWPT method. To make things simple we assume that the Gaussian plume is within a square two- dimensional domain and that the maximum concentration, C0, occurs at the center of the domain. In this case, the concentration as function of position and time is given by C(x, y, t) = C0s 2 0 s2 e −(x−x0)2 2s2 − (y−y0)2 2s2 (2.20) To calculate a reasonable value for the local-dispersion coefficient we assume a pore water velocity v = 10−7m/s and dispersivity α = 1 cm which results in a local-dispersion 42 10−3 10−2 10−1 10−7 10−6 10−5 10−4 10−3 10−2 10−1 h Er ro r MC σ/Δx = 1 WC σ/Δx = 1 MC σ/Δx = 0.3 WC σ/Δx = 0.3 MC σ/Δx = 0.1 WC σ/Δx = 0.1 σ/dx = 0Δx/h=1.0 10−3 10−2 10−1 10−7 10−6 10−5 10−4 10−3 10−2 10−1 h Er ro r MC σ/Δx = 1 WC σ/Δx = 1 MC σ/Δx = 0.3 WC σ/Δx = 0.3 MC σ/Δx = 0.1 WC σ/Δx = 0.1 σ/Δx = 0Δx/h=0.66 Figure 2.2: Error for one-dimensional simulations as a function of the smoothing length h for Monte Carlo (MC) and Weakly Compressible (WC) formula- tions defined by equations (2.8) and (2.9), respectively. Error magnitude is shown for different ratios of particle spacing over smoothing length ∆x/h and for different perturbations over particle spacing, σ/∆x. Uni- form particle spacing corresponds to σ/∆x = 0. 43 MC-SPH RWPT Particles carry solute concentration Particles carry solute mass It is possible to compute chemical reactions at individual particles Chemical reactions must be evaluated at some cell scale It is more complex to implement full dispersion tensor It is easy to implement full dispersion tensor Numerical precision to represent concentration values up to hardware representation Numerical precision to represent concentration values given by number of particles Simulates solute mass transfer between particles Simulates solute mass transfer at cell scale It demands more computational effort but it is possible to get very good accuracy with moderate use of memory It demands less computational effort but it requires more memory to get higher accuracy Table 2.1: Comparison of MC-SPH and RWPT methods. coefficient equal to D = v · α = 10−9m2/s. We consider that the plume is centered in a square domain of side L = 100m and that s0 = 5m, so that boundary effects are negligible. To integrate the solution in time we use a time step ∆t = 11.6 days and we consider a total simulation time T = 500∆t = 15.9 years. We performed a series of simulations to evaluate the relative performance and conver- gence properties of the RWPT and our MC-SPH methods. First, we solved the problem using a RWPT method with different combinations of averaging volumes and number of particles to represent the mass in the cell with the maximum concentration as summa- rized in Table 2.2. Then, we simulated the same situation using MC-SPH considering different combinations of number of particles and kernel smoothing length which defines the average number of particles per kernel support volume as summarized in Table 2.3. We note that the direct comparison of both methods is difficult because of the differences in the way they calculate concentrations that result in different spatial resolutions and accuracies for the same number of particles. For example, particles in the RWPT simu- lations are only located within the plume edges while in MC-SPH simulations particles, as explained below, must cover a larger area. In the MC-SPH simulations presented here particles are quasi-randomly distributed in all the domain. However, it would be possible to improve the spatial resolution by distributing the same number of particles in a smaller area. Nevertheless, we believe that the results presented next constitute a fair comparison of both methods. 44 Simulation # cells Nr Np RW1 50 x 50 10 329 RW2 50 x 50 100 3,957 RW3 50 x 50 1000 40,713 RW4 100 x 100 10 1,257 RW5 100 x 100 100 15,345 RW6 100 x 100 1000 157,972 RW7 200 x 200 10 4,993 RW8 200 x 200 100 60,965 RW9 200 x 200 1000 627,153 Table 2.2: Parameters used in RWPT simulations: number of cells used to calculate concentrations, number of particles used to represent the mass within the cell with maximum concentration Nr, and total number of particles Np. Simulation h Np Nk SPH1 0.5 10,000 3 SPH2 0.5 20,000 6 SPH3 0.5 40,000 13 SPH4 0.5 60,000 19 SPH5 0.5 80,000 25 SPH6 1.0 10,000 13 SPH7 1.0 20,000 25 SPH8 1.0 40,000 50 SPH9 1.0 60,000 75 SPH10 1.0 80,000 101 SPH11 2.0 10,000 50 SPH12 2.0 20,000 101 SPH13 2.0 40,000 201 SPH14 2.0 60,000 302 SPH15 2.0 80,000 402 Table 2.3: Parameters used in SPH simulations: smoothing length h, total number of particles Np, and average number of particles per kernel support volume Nk. 45 2.3.2.1 Initial particle and concentration distribution To apply the RWPT method to this problem we must first map the spatial concentration distribution to a regular Cartesian grid. Next, we must set the number of particles that represent a given mass to compute the equivalent particle distribution. In this example we use different numbers of particles to represent the mass contained in the cells with highest concentration value. Particles within each cell are initially distributed using a quasi-random distribution to generate a uniform spatial coverage. Figure 2.3 shows the corresponding particle distribution and the equivalent cell concentrations for some exam- ple configuration. We note that particles are only present in cells where concentration values are greater than some numerical threshold equal to the ratio between the mass of individual particles and the cell volume (see Appendix C for details). Because of the av- eraging procedure used to compute cell concentrations, the maximum cell concentration value is less than C0. Similar differences between cell values and actual concentrations occur in the rest of the domain and they are relatively more important near the plume edges where cells contain fewer particles. Figure 2.3: In RWPT simulations concentration values are approximated according to the spatial distribution of particles. The left figure show the initial particle distribution corresponding to a Gaussian plume with maximum concentration C0 at the domain center. The right figure shows concen- tration values computed according to the number of particles in each cell. The initialization of particle positions and concentration values in SPH simulations are independent. Particle positions are assigned such that the resulting particle distribution covers the region of interest. For example, given an initial Gaussian plume and the same number of particles as used in Figure 2.3, particles can be quasi-randomly distributed in a 46 rectangular volume or distributed in a uniform rectangular lattice within a circular region as shown in Figure 2.4. We observe that the maximum concentration value is within 0.1% of the actual maximum concentration value C0 in both cases. This shows that the MC- SPHmethod provides a better numerical resolution to represent concentration values than the RWPT using the same number of particles. In this simple example, particles are only created within a region of the numerical domain where concentrations are greater than a given threshold value plus a surrounding buffer zone. The buffer zone is necessary to provide additional space to allow dispersion to distribute the initial solute mass in a larger volume. In real simulations it is important to prevent the existence of isolated particles with non-negligible concentration at the edge of the particle cloud to avoid numerical errors. There are three possible alternatives to control this source of error: (i) generate a particle distribution that covers all the domain, (ii) use a dynamical scheme that inserts particles as needed, (iii) generate an initial particle distribution with a buffer zone large enough to guarantee an appropriate particle distribution during the simulation. The first alternative is very simple to implement but it can become prohibitively expensive in large-scale simulations or in simulations that require fine-spacing between particles. The second solution works well but it is more difficult to implement and introduces some computational overhead because it requires more sophisticated data structures to store and manage the particle set. The third alternative combines the advantages of the other two because it is very easy to implement and requires fewer particles than the first option. Figure 2.4: In MC-SPH simulations concentration values are directly assigned to each particle. Figures show two possible initial particle distributions and cor- responding concentration values equivalent to a Gaussian plume with maximum concentration C0 at the domain center. 47 2.3.2.2 Performance The overall number of floating operations and memory requirements of the RWPT method using a background grid to compute concentration scales linearly with the to- tal number of particles, i.e. it is O(N). However, the application of (C.3) to compute smoother concentration distributions requires O(N2) operations. On the other hand, the evaluation of the summation in (2.6) corresponds to an n-body problem which naive implementation scales as O(N2) (Greengard, 1994). However, because the compact sup- port of the kernel the actual work required to compute the summation can be reduced to O(NNk), where Nk corresponds to the average number of particles per kernel sup- port volume. The implementation of such algorithms requires an efficient method for searching near-neighbor particles. Such algorithms are based on data structures used to classify particles according to their spatial coordinates. For constant smoothing length implementations as presented in this paper, the background grid algorithm is the most efficient method (Viccione et al., 2008). For spatially varying smoothing lengths, more sophisticated data structures based on octrees or binary trees must be used (e.g. Barnes and Hut, 1986; Waltz et al., 2002). An explicit implementation of the proposed meshless method as discussed in Section 2.2.1 requires an amount of memory that scales with the number of particles (O(N)). Figure 2.5 shows the CPU time required to complete a single time step as function of the total number of particles in RWPT simulations. As expected the computational cost of the method grows linearly with the total number of particles. On the other hand, Figure 2.6 shows that in SPH simulations the CPU time depends in a non-linear fashion on the total number of particles and the kernel smoothing length. Larger smoothing lengths, equivalent to more particles per kernel support volume, result in longer simulation times for the same total number of particles. Figure 2.6 also shows that, as expected, the CPU time required to complete a single time step scales linearly with the product of the total number of particles and the average number of particles per kernel support, NPK = NpNk. The differences observed between the curves corresponding to different smoothing lengths for the same product NPK are due to differences in performance of the routine that evaluates the changes in concentration at each time step as result of different combinations of Np and Nk. 48 0 1 2 3 4 5 6 7 x 105 0 200 400 600 800 1000 1200 Np No rm ali ze d CP U Ti m e Figure 2.5: Normalized CPU time required to solve one time step using RWPT as function of the total number of particles Np. Computational cost of RWPT simulations is proportional to the total number of particles Np. 49 1 2 3 4 5 6 7 8 x 104 0 50 100 150 200 250 300 Np No rm ali ze d CP U Ti m e h = 0.5 h = 1.0 h = 2.0 0 0.5 1 1.5 2 2.5 3 3.5 x 107 0 50 100 150 200 250 300 Np Nk No rm ali ze d CP U Ti m e h = 0.5 h = 1.0 h = 2.0 Figure 2.6: Normalized CPU time required to solve a single time step as function of the total number of particles Np, kernel smoothing lenght h, and average number of particles per kernel support volume Nk. Computational cost of MC-SPH simulations is proportional to the product of Np and Nk. 50 2.3.2.3 Accuracy We used two criteria to compare the accuracy and convergence properties of RWPT and MC-SPH. Figure 2.7 compares the temporal evolution of the simulated maximum concentration with the theoretical result according to (2.20). The simulated RWPT results exhibit unphysical oscillations and large errors. Such unphysical oscillations would create serious problems if an operator splitting approach was used to simulate reactive transport simulations where these errors would be amplified by non-linear reactions. Local errors in RWPT simulations do not only depend upon the total number of particles but also on the number of particles at each cell and the cell volume. For example, simulation RW3 with Np = 40713 produces maximum concentration values that are closer to the true value than the results of RW6 and RW9 with Np=157972 and Np = 627153, respectively. Simulations with lower Nr as defined in Appendix C, not shown in the figure, produced results with even larger errors. Figure 2.7 also shows a comparison of the maximum concentration simulated with MC-SPH considering Np = 40, 000 and the analytical solution. All of the simulations give solutions that are free of unphysical oscillations, however, simulations with low average number of particles per kernel support volume such as SPH3 (Nk = 13) can result in considerable local errors. Local errors can be made negligible by choosing an appropriate combination of total number of particles and kernel smoothing length to obtain larger average number of particles that effectively contribute to the numerical integration, e.g. simulations SPH8 (Nk = 50) and SPH13 (Nk = 201). Results of other simulations with larger number of particles, not shown in the figure, produced smaller errors. Figure 2.8 shows the global error as function of the total number of particles in the simulation and CPU time per time step. It is clear that using these metrics numerical solutions computed using MC-SPH converge faster to the true solution than the RWPT solutions. Figure 2.8 also shows that MC-SPH solutions have different convergence rates depending upon the kernel smoothing length used. Moreover, curves corresponding to different smoothing lengths intersect indicating the transition between regions where the error is controlled by the average number of particles per kernel support volume Nk (small Np) and the region where the error depends upon the spatial resolution given by the kernel smoothing length (large Np). For simulations that require low CPU time, the convergence rate for MC-SPH simulations is faster than the one for RWPT simulations. However, the convergence rate of both methods becomes similar for simulations with larger number of particles that require longer CPU times to complete a single time step. 51 0 50 100 150 200 250 300 350 400 450 5000.6 0.7 0.8 0.9 1 Time Step C m ax /C 0 Analytical RW3 RW6 RW9 0 50 100 150 200 250 300 350 400 450 5000.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Time Step C m ax /C 0 Analytical SPH3 SPH8 SPH13 Figure 2.7: Maximum simulated concentration versus time step. RWPT simulations with resolution number Nr = 1000 and MC-SPH simulations with Np = 40000. Estimated concentrations in RWPT simulations present numerical oscillations due to representing the solute mass distribution as a finite set of particles. MC-SPH simulations compute concentrations that are free of numerical oscillations. 52 102 103 104 105 106 10−6 10−5 10−4 10−3 Np Er ro r RWPT SPH h = 0.5 SPH h = 1.0 SPH h = 2.0 100 101 102 103 104 105 10−6 10−5 10−4 10−3 Normalized CPU Time Er ro r RWPT SPH h= 0.5 SPH h = 1.0 SPH h = 2.0 Figure 2.8: Normalized global error versus total number of particles and normalized CPU time. Error computed as defined in (2.19) substituting N by the number of cells for RWPT and the total number of particles for MC-SPH simulations. The convergence rate, measured as function of the total number of particles Np or the CPU time required to complete a single time step, is faster for MC-SPH simulations than for RWPT simulations. However, the convergence rate of both methods become similar for simu- lations that use more particles and demands longer CPU times. 53 2.3.3 Advection–Dispersion in Heterogeneous Porous Media The objective of this section is to evaluate the performance of the MC-SPH approach to simulate solute transport in two-dimensional heterogeneous porous media. To verify our MC-SPH code we compared it with the ULTIMATE total variation diminishing (TVD) finite difference solver and the hybrid method of characteristics (HMOC) particle-mesh solver included in the popular MT3DMS package (Zheng and Wang, 1999). We focus our analysis on verifying if the numerical results satisfy some basic physical requirements such as: avoiding numerical dispersion, providing positive solution free of oscillations, and mass conservation. 2.3.3.1 Setup Before solving the solute transport problem we generated a velocity field as follows: (i) generate a moderately heterogeneous random hydraulic conductivity field, (ii) calculate a velocity field by solving the saturated flow problem using MODFLOW (Harbaugh, 2000) considering a constant hydraulic head gradient from left to right, and no-flow boundaries at top and bottom. Table 2.4 shows the parameters used to generate the random hydraulic conductivity field and to solve the flow problem. We used the resulting velocity field to simulate conservative transport of a square initial plume with constant concentration, C0. Table 2.5 shows the parameters used to solve the transport problem. Figure 2.9 shows a schematic of the simulation setup. In all the simulations discussed below we only considered constant isotropic dispersion. Description Symbol Value Variance of ln(K) σY 1.4 Correlation length of ln(k) IY 2.5 Domain dimension (Lx, Ly) (200IY , 50IY ) Grid size ∆ IY /5 Mean velocity in x U 0.81 Max. velocity in x umax 6.22 Table 2.4: Parameter and results of flow model. For the TVD simulations we used the same grid discretization that was used to solve the flow problem. For the HMOC simulations the allowed total maximum number of 54 Description Symbol Value Initial plume size (Lpx, Lpy) (20IY , 20IY ) Initial plume center (Xp, Yp) (36IY , 25IY ) Péclet number Pe = UIY /D [20, 200,∞] Dimensionless time step τ = U∆t/∆ 6.5 · 10−3 Mean CFL number CFLmean = U∆t/∆ 0.03 Max. CFL number CFLmax = umax∆t/∆ 0.25 Table 2.5: Parameter values used in transport model. Figure 2.9: Domain dimensions, square initial plume, and breakthrough observation points P1 and P2 along centerline. 55 particles was 6 · 105 and the number of particles per cell in cells where the relative concentration gradient (Zheng and Wang, 1999, p. 65) was higher than 1 · 10−5 was 15. The threshold value, DCHMOC, which controls if the forward or backward MOC method is used for an individual cell according to its relative concentration gradient was set equal to 1 · 10−4 (Zheng and Wang, 1999, p. 73). In the MC-SPH simulations the total number of particles was constant during the simulation and equal to 3.8 · 105 which was equivalent to an initial number of particles per grid cell equal to 8. Particles where initially distributed in a rectangular lattice within a rectangular region of size 68IY in the direction of the flow and 28IY in the transverse direction centered at the initial plume center. The kernel smoothing length was constant and equal to half the grid size used in the MT3DMS models, so the spatial resolution of the three methods was comparable. The results of both, TVD and HMOC, methods correspond to concentration values at the center of grid cells. We interpolated the MC-SPH results which correspond to con- centration values at scattered points onto a similar grid to compare them. We stress that this interpolation was needed only for comparison purposes and it is usually not necessary in SPH simulations. We used the following expression to compute interpolated values AI , AI (ri) = j AjŴ (|ri − rj|) j Ŵ (|ri − rj|) = j Ajψij (2.21) where ψij are Shepard functions (Shepard, 1968) and summations are over all particles. Although the interpolation in (2.21) is valid for any kernel function Ŵ , we used the same kernel used in MC-SPH simulations to get interpolated quantities with similar spatial resolution. 2.3.3.2 Results Figure 2.10 shows the spatial concentration distribution at the end of the simulation for the advection-only case. Solutions given by the three methods are very different. The TVD solver is not able to avoid numerical dispersion so the plume exhibits large dilu- tion and the initial plume mass is distributed within a much larger volume. There are only few areas where the solute concentration is similar to the initial concentration. The 56 HMOC produces less numerical dispersion and the plume edges are clearly distinguish- able. However, the effects of numerical dispersion are clear in zones located between fast or slower fingers and in front and behind the main plume. The existence of those arti- ficial low concentration zones has some important practical implications. For example, in presence of reactions controlled by mixing such as biodegradation those low concen- tration zones could potentially extend the reactive zone near the plume edge (Cirpka et al., 1999b). The concentration distribution generated by the MC-SPH code is free of numerical dispersion and the plume edges are very sharp as expected in absence of local-dispersion. Zones without contaminant located between fast and slow fingers are observable surrounding all the plume volume. It is interesting to notice some isolated high concentration spots in the front and back of the plume as result of connected high and low permeability regions. Figure 2.11 shows concentration values along the domain centerline at the end of the simulation, i.e., after the plume center has traveled 62 integral scales. The three methods produce very similar results for low Péclet numbers. As expected, the results produced by TVD and HMOC methods are identical considering that both methods share the same dispersion solver routine. It is more interesting to notice the good agreement between MC-SPH and the other two methods for low Péclet values. It is difficult to say if the small differences observed at the plume edges are due to differences between the meth- ods or the interpolation method used to map the MC-SPH results onto a grid. More important differences are observed for the more strongly advection-dominated scenarios. For Péclet number equal to 200, the two Lagrangian based methods, HMOC and MC- SPH, perform similarly while the TVD solution is smoothed by numerical dispersion. This comparison clearly shows that even high-order Eulerian mesh-based methods such as TVD cannot compete with particle-based methods for advection-dominated problems. For the advection-only case the three methods give solutions that are clearly distinguish- able. TVD results show little difference with respect to the situation for Pe = 200. On the other hand, HMOC and MC-SPH results are closer to the expected sharp profile with concentration values equal to the initial concentration or zero. The MC-SPH solu- tion seems to perform better close to the plume edges where the HMOC solution gives a smoother profile as consequence of the accumulated numerical dispersion due to the interpolation of concentration values from scattered points to the cell centers at each time step. Figures 2.12 and 2.13 show breakthrough curves for points P1 and P2 located along the domain centerline at 26IY and 42IY downstream from the initial plume center, respec- 57 Figure 2.10: Spatial concentration distribution at dimensionless time τ = Ut/IY = 62 for Pe = ∞. TVD solver (top), HMOC solver (middle), and SPH solution mapped onto rectangular grid (bottom). Only dimensionless values C/C0 > 0.001 are shown. 58 0 10 20 30 40 50 60 70 800 0.2 0.4 0.6 0.8 1 C / C 0 Pe = 20 0 10 20 30 40 50 60 70 800 0.2 0.4 0.6 0.8 1 C / C 0 Pe = 200 0 10 20 30 40 50 60 70 800 0.2 0.4 0.6 0.8 1 d / IY C / C 0 Pe = ∞ TVD HMOC SPH Figure 2.11: Concentration versus accumulated distance along centerline at dimen- sionless time τ = Ut/IY = 62. 59 tively. Breakthrough curves for low Péclet values equal to 20 given by the three methods are almost identical. This confirms that the MC-SPH method is able to simulate sit- uations where dispersion is important with accuracy that is comparable to well-tested mesh-based methods as used in MT3DMS. As observed in the longitudinal profile com- parison, the difference between HMOC and MC-SPH solutions and the TVD solution increases as the transport process becomes controlled by advection. For Péclet number equal to 200, HMOC and MC-SPH produces similar results while the corresponding TVD breakthrough curves have a completely different shape typical of much higher dispersion coefficients. Finally, for the purely-advective case the breakthrough curves corresponding to each method become very different. For the point located closer to the plume center, P1, HMOC and TVD predict a much earlier arrival time and longer tail. The earlier arrival time is due to the lateral mixing produced by the numerical dispersion, which transfers solute concentration from faster plume fingers that pass close to point P1. The longer tail is also due to lateral and longitudinal numerical mixing. The breakthrough curve corresponding to the MC-SPH solution is not affected by numerical dispersion, thus it exhibits a rectangular shape as expected. For point P2, the three breakthrough curves have similar arrival time indicating that a set of fast streamlines forming a fast advancing front passes through P2, so the lateral mixing does not have a significant ef- fect. However, lateral and longitudinal num# # There is insufficient memory for the Java Runtime Environment to continue. # Native memory allocation (mmap) failed to map 135266304 bytes for committing reserved memory. # An error report file with more information is saved as: # /usr/local/web-apis/harvester/bin/hs_err_pid61754.log
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Particle and streamline numerical methods for conservative...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Particle and streamline numerical methods for conservative and reactive transport simulations in porous… Herrera, Paulo Andres Ricci 2009
pdf
Page Metadata
Item Metadata
Title | Particle and streamline numerical methods for conservative and reactive transport simulations in porous media |
Creator |
Herrera, Paulo Andres Ricci |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | Reactive transport modeling has become an important tool to study and understand the transport and fate of solutes in the subsurface. However, the accurate simulation of reactive transport represents a formidable challenge because of the characteristics of flow, transport and chemical reactions that govern the migration of solutes in geological formations. In particular, solute transport in natural porous media is advection-controlled and dispersion is higher in the direction of flow than in the transverse direction. Both characteristics create difficulties for traditional numerical schemes that result in numerical dispersion and/or spurious oscillations. While these errors can often be tolerated in conservative transport simulations, they can be amplified in presence of chemical reactions resulting in much larger errors or unstable solutions. In this thesis, new Lagrangian based methods to simulate conservative and reactive transport in porous media are investigated. First, the derivation of a new meshless approximation based on smoothed particle hydrodynamics (SPH) to simulate conservative multidimensional solute transport, including advection and anisotropic dispersion, is presented. Second, a hybrid scheme that combines some of the advantages of streamline-based simulations and meshless methods and that allows simulating longitudinal and transverse dispersion without requiring a background grid is also derived. The numerical properties of both methods are analyzed analytical and numerically. Furthermore, both formulations are compared with existing numerical techniques in a set of two- and three-dimensional benchmark problems. It is demonstrated that the proposed schemes provide accurate and efficient solutions of physical transport processes in heterogeneous porous media and overcome most of the issues in existing numerical formulations. The new methods have the potential to remove or minimize numerical dispersion and grid orientation effects and, in the case of the hybrid streamline method, also eliminate spurious oscillations even in presence of large longitudinal to transverse dispersivity ratios. Therefore, the results presented in this thesis confirm that the Lagrangian formulations of solute transport investigated here are viable and compelling alternatives to simulate reactive transport versus more standard numerical techniques. |
Extent | 10878280 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-11-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0053008 |
URI | http://hdl.handle.net/2429/15967 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geological Engineering |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_spring_herrera_paulo.pdf [ 10.37MB ]
- Metadata
- JSON: 24-1.0053008.json
- JSON-LD: 24-1.0053008-ld.json
- RDF/XML (Pretty): 24-1.0053008-rdf.xml
- RDF/JSON: 24-1.0053008-rdf.json
- Turtle: 24-1.0053008-turtle.txt
- N-Triples: 24-1.0053008-rdf-ntriples.txt
- Original Record: 24-1.0053008-source.json
- Full Text
- 24-1.0053008-fulltext.txt
- Citation
- 24-1.0053008.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0053008/manifest