Numerical Modeling Techniques for Assessing Three-Dimensional Diffusion Processes in Heterogeneous Rock Samples on the Sub-mm Scale by Scotland Russell Ellis B.A.Sc., The University of British Columbia – Vancouver, 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Geological Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2010 © Scotland Russell Ellis, 2010 ABSTRACT The burial of high level nuclear wastes in geologic repositories requires careful consideration of the long-term hydrogeological and geochemical stability of the receiving formations. Care must be taken due to the high environmental sensitivity of the waste material and its potential long-term effect on groundwater. Studies into the host rocks' natural ability to minimize contaminant migration are a matter of high priority in the planning for long-term storage of high level radioactive waste. Focusing on porous sedimentary rock, this study aims to examine numerical techniques used to analyze and interpret experimental data that characterize the distribution of porosity in geologic samples at the µm to mm scale. Because repositories are hosted in natural fine-grained rock formations, transport in the vicinity of these repositories is diffusion-controlled and believed to be affected substantially by heterogeneities at these scales. Data available for the analyses consist of non-destructive, high-resolution measurements of porosity obtained using new developments in X-ray imaging. Advances in computing technology make it possible to numerically analyze the expected patterns of sub- mm-scale diffusive transport for these large experimental data sets. The modeling analyses examine 3D diffusive transport in heterogeneous rock samples and evaluate the effect of data resolution and image processing techniques on the connectivity of the transport pathways. The simulation results provide insight into small- scale diffusive transport of solutes, and guide the needs for dataset resolution and handling for these types of investigations. With increased availability of experimental results, further modeling studies could be conducted. These studies would aim at developing a link between simulation results and observed data to further develop the transport theory for contaminant migration on this scale. ii TABLE OF CONTENTS ABSTRACT ............................................................................................................................ ii TABLE OF CONTENTS ...................................................................................................... iii LIST OF TABLES ..................................................................................................................vi LIST OF FIGURES .............................................................................................................. vii ACKNOWLEDGEMENTS ....................................................................................................x Chapter 1: Introduction ..........................................................................................................1 1.1 Background ...................................................................................................................1 1.2 Transport Parameters and Processes .........................................................................2 1.3 Previous Work ..............................................................................................................7 1.4 Motivation for Current Project .................................................................................10 1.5 Thesis Objectives ........................................................................................................11 Chapter 2: Methods ...............................................................................................................13 2.1 Introduction.................................................................................................................13 2.2 Data Collection and Format ......................................................................................13 2.2.1 Data Collection ......................................................................................................13 2.2.2 Data Cropping........................................................................................................15 2.2.3 Data Binning..........................................................................................................16 2.3 Modeling Methods ......................................................................................................17 2.3.1 Modeling Domain .................................................................................................17 2.3.2 Boundary and Initial Conditions .........................................................................17 2.3.3 Diffusion Modeling...............................................................................................17 2.3.4 Inter-Facial Fluxes - Averaging Formulations...................................................18 iii Chapter 3: MIN3P Diffusion Model Verification ...............................................................22 3.1 Introduction..................................................................................................................22 3.2 Verification ...................................................................................................................22 3.2.1 Diffusion in 1D Semi-Infnite Domain .................................................................23 3.2.2 Diffusion in 2D Semi-Infinite Domain with Patch Source.................................25 3.2.3 Diffusion in 3D Infinite Domain with Interior Brick Source.............................29 3.3 Diffusion Number Analyses ........................................................................................31 3.3.1 1D Time Step Analysis..........................................................................................32 3.3.2 3D Time Step Analysis..........................................................................................34 Chapter 4: Sub-mm Scale Diffusion Modeling in a Sandstone Core ................................38 4.1 Introduction..................................................................................................................38 4.2 Dataset Description......................................................................................................39 4.2.1 Cropping and Binning ...........................................................................................40 4.2.2 Binary Datasets ......................................................................................................42 4.2.3 16 Bit (Porosity Map) Datasets..............................................................................44 4.3 Diffusion modeling.......................................................................................................45 4.3.1 Input Parameters .................................................................................................46 4.4 Results ...........................................................................................................................48 4.4.1 Simulations of 3D Diffusion for Binary Dataset.................................................48 4.4.2 Crop 1 ....................................................................................................................54 4.4.2.1 Binary Dataset ................................................................................................54 4.4.2.2 16 Bit Dataset .................................................................................................56 4.4.2.3 Binary Versus 16 Bit Datasets........................................................................60 4.4.2.4 Arithmetic Versus Harmonic Averaging........................................................61 iv 4.4.3 Crop 2 ....................................................................................................................63 4.4.3.1 Binary Dataset ................................................................................................63 4.4.3.2 16 Bit Dataset .................................................................................................66 4.4.3.3 Binary Versus 16 Bit.......................................................................................66 4.4.4 Crop 1 Versus Crop 2 ...........................................................................................68 4.5 Discussion ....................................................................................................................69 Chapter 5: Summary and Conclusions................................................................................73 References...............................................................................................................................78 Appendices .............................................................................................................................81 v LIST OF TABLES Table 3.1. Simulation parameters used for 1D diffusion verification ....................................24 Table 3.2. Simulation parameters used for 2D diffusion verification. ...................................26 Table 3.3. Simulation parameters used for 3D diffusion verification. ...................................29 Table 3.4. Simulation parameters used for 1D diffusion verification ....................................32 Table 3.5. Simulation parameters used for 1D diffusion verification ....................................35 Table 4.1. Summary of porosity data sets created for simulations.........................................39 Table 4.2. Binning of binary dataset crop 1............................................................................43 Table 4.3. Binning of 16 bit dataset crop 1.............................................................................45 Table 4.4. Major input parameters for sandstone core diffusion simulations ........................47 Table 4.5. Binned parameters for binary dataset crop 1 .........................................................54 Table 4.6. Binned parameters for 16 bit dataset crop 1 ..........................................................57 Table 4.7. Comparison of ratios of porosity / mass in............................................................60 Table 4.8. Parameters of arithmetic vs harmonic ...................................................................61 Table 4.9. Parameters for binary dataset crop 2 as a function of binning degree...................64 Table 4.10. Binned parameter for 16 bit dataset crop 2..........................................................66 Table 4.11. Comparisons between crops for unbinned datasets.............................................68 vi LIST OF FIGURES Figure 1.1. Below the REV scale, randomly selected small volumes V1, do not characterize the entire system. Above the REV scale, V2 has a porosity that is representative of the material ......................................................................................................................................4 Figure 1.2. Conceptual model of the tortuosity of a material. Tortuosity is calculated using the differences between a straight line path (L) and the actual effective path taken by a particle (Le) ................................................................................................................................6 Figure 2.1. Size and orientation of core sample .....................................................................14 Figure 2.2. a) 2D image slice of 11mm diameter sandstone taken from 3D dataset obtained using MicroCT scanner. b) Section of sample showing binary spatial distribution of porosity where the white regions correspond to open pores. Ratio of white pixels to total pixels corresponds to average porosity of sample..............................................................................15 Figure 2.3. a) Binary porosity dataset with 144x144 voxels. Binary 8 bit histogram shows values at 0 and 255 b) Same porosity dataset after being binned 2x (to 72x72 voxels). The histogram now shows a spread in porosity values...................................................................16 Figure 2.4. a) Arithmetic averaging will overestimate fluxes into areas of low porosity, b) Harmonic averaging between the cells ensures that large fluxes are not present between cells of high and low porosity ..........................................................................................................19 Figure 3.1. Diffusion at multiple output times along 1D column. MIN3P simulation (points) and Ogata-Banks analytical solution (solid lines) ...................................................................25 Figure 3.2. a) 2D patch source model setup. Patch source is located at x=0 and is infinite in z with a finite width in the y-direction of 7.5% the extent of the y domain. b) Diffusion contours in the x-y plane at T = 50 days..................................................................................27 Figure 3.3. Diffusion after 50 Days with multiple diffusion coefficients from a 2D patch source. MIN3P numerical simulations (points), Modflow numerical simulations (dashed lines) and Patch3D analytical solution (solid lines) ................................................................28 Figure 3.4. 3D brick source model setup. The source is located in the middle of the model domain and is in the shape of a brick with dimensions (x,y,z) of 0.009,0.003,0.003 metres..29 Figure 3.5. The brick provides a source of finite solute mass from which diffusion occurs in all three dimensions in an expanding oblate spheroid, concentration contours shown at T = 50 days ..........................................................................................................................................30 Figure 3.6. 3D diffusion outward from a brick source using multiple diffusion coefficients. MIN3P numerical simulations (symbols) and brick analytical solution (solid lines), T = 50 days ..........................................................................................................................................31 vii Figure 3.7. Percent error after 50 days along 1D column. MIN3P simulations using differing maximum time steps. ...............................................................................................................33 Figure 3.8. Concentration profile after 50 days along 1D column. MIN3P simulations using differing maximum time steps. ................................................................................................34 Figure 3.9. Percent error after 0.1 Hours along a 3D column. MIN3P simulations using differing maximum time step...................................................................................................36 Figure 3.10. Contour plot of concentrations after 0.1 hours using time steps of 0.0001 (black lines) and 0.01 hours (white lines). Cross section is through the middle of the column with source located at z = 0 m. ........................................................................................................36 Figure 4.1. Two separate 2x106 voxel datasets cropped from the original 1x108 voxel MicroCT porosity dataset. The top and bottom binary datasets have overall average porosities of 10% and 22% respectively. .................................................................................................41 Figure 4.2. Image of porosity of first slice and averaged porosities per slice versus distance along column. a) Crop 1 of the binary dataset. Average porosity of 22% b) Crop 2 of the binary dataset. Average porosity of 10%................................................................................43 Figure 4.3. Image of porosity of first slice and averaged porosities per slice versus distance along column. a) Crop 1 of the 16 bit dataset. Average porosity of 16.5% b) Crop 2 of the 16 bit dataset. Average porosity of 13.5%..............................................................................44 Figure 4.4. Model domain setup for numerical through diffusion experiment with Iodide reservoir input at z = 0. ............................................................................................................46 Figure 4.5. Concentration map at xy slice 19 from the source (z=3.3x10-4m) after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the binary dataset. ..............................49 Figure 4.6. Concentration map at xz slice along the centerline of y after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the binary dataset. ..............................................50 Figure 4.7. Concentration map after 1hr through xy slices a) 13 b) 26 c) 39 d) 52 e) 65 f) 78 and g) 91 out of 96 for the binary dataset................................................................................51 Figure 4.8. Concentration map at xy slice 19 from the source (z=3.3x10-4m) after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the dataset that has been binned by 4x. .................................................................................................................................................52 Figure 4.9. Concentration map at xz slice along the centerline of y after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the dataset that has been binned 4x...................53 Figure 4.10. Concentration map at an xy slice 3.3x10-4m from the source after 0.01hrs for the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x.................................................55 viii Figure 4.11. Concentration profiles along the long axis of the domain for each output time for the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x ..........................................56 Figure 4.12. Concentration map at an xy slice 3.3x10-4m from the source after 0.01hrs for the porosity dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x..............................................58 Figure 4.13. Concentration profiles along the long axis of the domain for each output time for the porosity map dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x ...............................59 Figure 4.14. Slice of porosity 3.3x10-4 m from the source for binary data set, as well as concentration maps for models with harmonic and arithmetically averaged modeling schemes (respectively) ...........................................................................................................................62 Figure 4.15. Slice of porosity for 16 bit data set (porosity map), as well as concentration maps for models with harmonic and arithmetically averaged modeling schemes (respectively) .................................................................................................................................................62 Figure 4.16. Concentration map for crop 2 at an xy slice 3.3x10-4m from the source after 0.01hrs for crop 2 of the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x...............64 Figure 4.17. Concentration profiles for crop 2 along the long axis of the domain for each output time for the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x.......................65 ix x ACKNOWLEDGEMENTS I would like to acknowledge all those who made my thesis work possible. Thanks to my advisor, Dr. Uli Mayer. Uli was a very patient and supportive guide and mentor through the years and it was privileged to work with him. Thanks to professors Dr. Roger Beckie and Dr. Leslie Smith for providing me with the bulk of my studies to keep me on my toes and for being there to help when it was needed. Thanks to my colleagues across the country at UNB, Dr. Tom Al and Mosi Agbogun for their contributions to the project that made my work possible. I would also like to thank my colleagues at UBC, particularly Katie Jones and Sergi Molins for their continual contributions and patience throughout the years. A very special thanks to my family for supporting me throughout the whole journey. This project was made possible by funding through an NSERC (Natural Sciences and Engineering Research Council of Canada) strategic project grant. Chapter 1: Introduction 1.1 Background Transport and containment of contaminants in the subsurface are issues that are of importance worldwide. Within our technological society, increasing energy demand creates a need for reliable energy production. Accounting for ~15% of the world’s electricity supply (IEA, 2007), nuclear power has been and will remain an important global energy source. These power plants continually create high-level nuclear waste (used fuel rods) and the long- term contamination potential associated with these radioactive materials creates a technical challenge and societal concerns. As a result, investigations into deep geologic repositories for the storage of nuclear waste are numerous and ongoing. The difficulties with the storage of high level wastes tend to stem from the extremely slow radioactive decay timelines involved. Secure and safe containment for upwards of 100,000 years is commonly envisioned requiring the design of multiple safety barriers to warrant long-term waste isolation. These barriers may include emplacement of the waste in closed off, multiple-walled containers. These containers are built to withstand corrosive and mechanical stresses found at depth within the earth (Finnish Energy Industries, 2007). In the case of deep geologic repositories, an additional level of safety is provided by emplacing the waste storage facility deep underground in hydrogeologically isolated low permeability units. For the unlikely case that all other safety barriers should fail, geologic media provide additional protection for radionuclide migration away from the repositories. Because of their potential for minimizing contaminant migration, low permeability geologic media are often considered ideal sites for waste containment. For example, in Canada, low permeability argillaceous limestones are currently under investigation as potential host rock for a deep geologic repository for low and intermediate nuclear waste 1 (Intera, 2006). With permeabilities on the order of 10-12-10-13 m/s (Yllera 2004), the dominant transport process in these formations is expected to be diffusion. With diffusion- dominated transport and negligible potential for advective transport, these sites are likely suitable for containment of high level wastes for time periods exceeding 100,000 years. For this containment to be reliable, the time scale and spatial distribution of diffusive transport must be thoroughly examined. In these low permeability media, the diffusion rate is thought to depend on the materials’ physical properties such as porosity and tortuosity. Given their potential for small grain sizes and small values of porosity, the potential for large tortuosities is prevalent in low porosity materials such as sedimentary rocks (Koponen et al., 1996). 1.2 Transport Parameters and Processes Physical properties of porous media have a great effect on the transport processes occurring within the materials. In order to better understand the movements of contaminants, it is important that knowledge of the underlying transport processes be obtained. This is a common approach to problems in groundwater investigations where advection-dispersion is the dominant transport process (Gelhar et al., 1992, Harvey and Gorelick, 2000). For these studies, porosities and permeabilities are examined and average parameters defined at the scale of interest are used to quantify transport processes. Diffusion-dominated solute transport in groundwater has received somewhat less attention, in particular the effect of sub- mm scale heterogeneities on solute migration and contaminant distribution has not been investigated in detail. If transport is fully diffusion-dominated, the complexities of diffusion may be over-simplified, and a better understanding may be needed. This is the type of transport regime that is typically prevalent in host rocks suitable for the storage and containment of nuclear wastes. 2 To help further examine the transport processes and their parameters, the scale at which they can be linked must be determined. Many transport processes can be examined at the pore scale. At this scale, distinction can be made between the fluid filled pores and the samples’ solid matrix. The solid matrix makes up the load-bearing structure of the sample while the fluid filled pores are the regions where the transport can occur. Due to difficulties in determining the small scale matrix and pore structures, larger scale links between averaged parameters and transport processes are examined through the use of a representative elementary volume. This scale dependence is inherent to the concept of a representative elementary volume (REV). The concept of an REV is that within a certain minimum volume of material, the characteristics of that elementary volume become representative of the characteristics of the volumes around it (Figure 1.1). This minimum volume is the REV. When examining transport processes such as diffusion that usually occur at very small scales, the parameters that describe it can be very susceptible to the heterogeneity in the material properties that are found at that scale. If the investigation occurs below the REV scale, these spatially varying parameters may or may not link to the accepted views of diffusion that our models are based upon. 3 100% V1 V2 Volume Figure 1.1. Below the REV scale, randomly selected small volumes V1, do not characterize the entire system. Above the REV scale, V2 has a porosity that is representative of the material In the presently widely accepted views, diffusion is based on the second law of thermodynamics and random Brownian motion and is generally modeled through Fick’s Second law 2 2 x cD t c (1.1) 4 where c is a solute concentration, t is time, D is the diffusion coefficient, and x describes a spatial coordinate. Movement of a solute due to diffusion is thus modeled to be based upon the concentration gradient that is created between the source of the solute and the system into which it is introduced. This movement is hindered to a degree determined by the deposition and subsequent rearrangement of the grains and pores of the material through which the solute is traveling. Lithification and cementation can also affect the pore structure and the connectivity of pores in a rock formation. The pore structure creates a microscopic obstacle course for the solute through which a straight line path is generally impossible. The amount of deviation that an average solute particle takes from the straight line path is referred to as the tortuosity of the material, and it is ultimately described through the use of pore water or effective diffusion coefficients. Being dependent on microscopic transport paths of connected porosity, determining the tortuosity of a material directly is not a simple feat. As such, tortuosity is usually determined empirically through the relationship )( 0DDD peff (1.2) (Boving and Grathwohl 2001, Nakashima 2000). Where Deff is the effective diffusion coefficient, is porosity, Dp is the pore water diffusion coefficient, and D0 is the molecular diffusion coefficient in water. Thus tortuosity τ (>1) is D0/Dp or the ratio of the diffusivity in free water to the bulk diffusivity in the porous medium. A tortuosity value greater than 1 corresponds to a ratio of solute path lengths within a material where the longer effective path length (Le) is the divided by the straight line path length (L). 5 L Le Figure 1.2. Conceptual model of the tortuosity of a material. Tortuosity is calculated using the differences between a straight line path (L) and the actual effective path taken by a particle (Le) Some authors prefer to define tortuosity as Dp/D0 (Bear, 1979, Schwartz and Zhang, 2003), but as Bear (1979) explains, this is merely a matter of definition [(L/Le)2 vs. (Le/L)2]. In fully saturated systems, τ has been estimated as (e.g.: Schwartz and Zhang. 2003). 3 11 (1.3) When combining these two definitions, a relationship between Deff and Do can be derived that exclusively depends on porosity. 6 3 4 3 1 0 )( D Deff (1.4) This relationship can be generalized by introducing a variable exponent m: meff D D 0 (1.5) This relationship is commonly known as Archie’s law (e.g. Boving and Grathwohl, 2001) and allows the expression of effective diffusion coefficients as a function of porosity and the exponent m, which can be considered a fitting parameter. While these laws have been used extensively to explain the transport due to diffusion on large scales far above the lower limit of REV applicability, diffusion on the generally smaller scales of fully diffusion-dominated systems is not as thoroughly investigated. Using Fick’s and Archie’s laws as starting points, further investigations into the processes of diffusion can be undertaken by looking at the material parameters that may influence the transport of contaminants in more detail than has been accomplished before. Starting at the sub-mm scale, a more fundamental understanding of the processes may be obtained which can be upscaled to function with larger systems. 1.3 Previous Work Advances in material characterization techniques have led to a better understanding of the sub-mm scale characteristics of porous material. Since the mid 1990s, X-ray Computed Tomography (CT) and Magnetic Resonance Imaging (MRI) have been used to obtain high resolution images of porous media. While the earlier studies focused on larger scale flow and transport processes (Greiner et al. 1997, Oswald et al. 1997), smaller scale studies of porosity, pore structure, and permeability were soon to follow (Ruiz de Argondona et al. 7 1999, Van Geet et al. 2003, Al-Raoush and Willson 2005). Ruiz de Argondona et al. (1999) were able to use X-ray CT imaging as a method to non-destructively determine the porosity profiles and other internal features of dolomite building stones. They found that they were able to examine the evolution of pores and the creation of a non-uniform porosity distribution by freeze thaw weathering. At an even smaller scale, Al-Raoush and Willson (2005) used synchrotron x-ray microtomography to examine the pore structures of porous media systems fabricated of glass beads with a minimum diameter of 0.123 mm. With the high resolution (<10 μm) imaging capabilities of the synchrotron, they were able to create mechanistic idealizations of the pore network and determine the scale at which the REV might apply. However, given the small size of clay particles (<2 μm), even this high resolution imaging technique is not able to determine the pore scale structure of clays and shales. With the arrival of high resolution images of porous media, diffusive transport processes have been examined with more care and at a higher resolution. Nakashima (2000, 2003) has used imaging techniques to obtain 3D representations of solute and porosity distributions enabling the determination of effective diffusion coefficients through numerical analysis. From high resolution images of a natural sandstone, Nakashima et al. (2004) were able to analyse the pore structure and diffusive pathways through a random walk computer simulation. Theoretical effective diffusion coefficients coincided well with experimentally derived values. Unfortunately, the method used requires an explicit determination of the pore network structure based on sub-REV scale data from relatively large grained media (0.1 mm) and thus is not directly applicable to the finer grained geologic media that is of greater relevance for diffusion-dominated transport. 8 In the course of examining the diffusive transport processes, numerous authors have attempted to find and test correlations between image-obtained values of porosity and the corresponding diffusion coefficients calculated through modeling. Using a simple diffusion model with a single diffusion coefficient determined from flow-through diffusion tests, Boving and Grathwohl (2001) found a broad correlation based on Archie’s law between porosity and tortuosity. Indeed, their work confirmed that the dominating parameter controlling effective diffusion is porosity. The broadness of the correlation could be caused by simplicity of the diffusion model as the correlation seemed to fit best for certain samples used in the experimental procedures. To further constrain the relationship between porosity and diffusion, several authors have attempted to determine the diffusion coefficients more accurately by gathering concentration profile data (Tidwell 2000, Polak et al. 2003, Altman et al. 2004, Wersin 2004). Again the relationships were found to be broad and ill-defined due in part to the simple diffusion models that were applied (Tidwell et al. 2000, Altman et al. 2004). Using a single-rate diffusion model in which the same diffusion coefficient was applied throughout the system, Tidwell et al. (2000) qualitatively examined the variations of the measured diffusion profile from the expected single-rate diffusion model. From the observed variations, they surmised that diffusion not only depends on the magnitude of porosity but also its spatial distribution. Along with a simple single-rate diffusion model, Altman et al. (2004) also examined the potential of a multirate diffusion model based on a lognormal diffusion coefficient distribution. Unfortunately, this model, while being less conservative, gave similar results to the single-rate model and was thus not advantageous. With these models, they also found a positive correlation between porosity and diffusion, but were unable to quantify it. 9 Observing complex 3D diffusion behaviour, they saw the need for a more thorough diffusion model that would be able to account for 3D spatially variable diffusion coefficients. High resolution imaging of porous media can potentially offer insight into the spatial movements of differing particles in porous media. Depending on the scale of the pores relative to the resolution of the imaging device, binary pore scale images (distinguishing pores from matrix) or non-binary averaged porosity maps can be created to signify the porous media’s pore structure. Also available is the ability to threshold a non-binary porosity map into high and low values, effectively creating a binary pseudo pore-scale image. It is hypothesized that this thresholding technique may properly simulate the pore scale structure of the media. 1.4 Motivation for Current Project Considering that there is a need to investigate the role of sub-mm scale heterogeneities in rock formations, adequate imaging techniques to characterize and represent the porosity distributions must be developed. In addition, if solute migration on this scale is to be analyzed using numerical models, suitable numerical formulations to describe diffusion in heterogeneous systems must be available. This study will focus on investigating the effect of using different representations of observed porosity (e.g, binary and porosity map representations) and the resolution of the data on model results. In addition, different numerical formulations for calculating the diffusive fluxes will be evaluated. This method evaluation will provide the foundation for future integration of 3D porosity and transient concentration data in rock samples. A future goal is to use time- dependent measurements of solute concentration distributions to determine spatially distributed diffusion coefficients via transient modeling simulations and inverse methods. In 10 addition, it can be envisioned that the role of chemical reactions within the porous media and feedback on porosity distributions can be assessed using non-destructive imaging techniques combined with numerical analysis. However, these types of advances will only be possible if reliable and representative imaging techniques and numerical representations of porosity distribution and diffusive fluxes will be available. 1.5 Thesis Objectives The objectives of this work include the use of sub-mm scale diffusive transport modeling based on observed porosity data to examine the sensitivities of model results to changes in data resolution as well as model formulation. Specifically, the transport model will be used to answer the following questions: How do the diffusion profiles differ between a binary pore-scale model and a REV scale model, using a REV scale transport model? How does the resolution of the dataset affect the modeling of small scale diffusive fluxes? Can MicroCT imaging techniques adequately capture the pore scale intricacies of a sandstone sample? How does the pseudo pore-scale model compare to the porosity map derived through investigative techniques? Is there a scale below which REV scale modeling appears to no longer adequately describe the transport of the contaminant? 11 What is the effect of using harmonic versus arithmetic spatial weighting of diffusion coefficients on the model results? To address these objectives, transport modeling was conducted based on high resolution porosity images of sandstone samples. The overarching objective of this work is to evaluate the current status and the suitability of imaging and numerical techniques to analyze spatially heterogeneous diffusive transport on the sub-mm scale. Regrettably, experimentally derived solute concentrations for quantitative comparisons were not available at the time of writing and thus these works are composed of a thorough qualitative analysis. 12 Chapter 2: Methods 2.1 Introduction In order to simulate microscopic transport in geologic media, representative model parameters must be obtained at the scale of investigation. This study aims to examine diffusion processes on the sub-mm scale through a network of pores that permeate the media and provide pathways through which contaminants can diffuse. Depending on the grain size of the media, these pores can be on a scale as small as of 10-9 meters. This type of resolution cannot be provided using available 3D imaging techniques and could not be resolved using numerical models. Data is therefore collected on a slightly larger scale and REV scale models are used to simulate the diffusion processes For this study, a data set from a clean sandstone was chosen as the primary sample of interest due to its large pore size (approximately 10-5 meters). 2.2 Data Collection and Format 2.2.1 Data Collection Core samples that were considered to be of interest to the study were collected by a third party and sent to the University of New Brunswick (UNB). At UNB, the geologic samples were examined using micro CT scanning equipment to obtain 3D images of porosity with voxel dimensions (x,y,z) on the order of 18 μm per side. For comparison, very fine sand grains have diameters in the range of 62.5 – 125 μm. Using methods developed by Glover et al. (2006, 2009), bead packs of this grain diameter may have effective pore radii of approximately 15 μm. Depending on cementation and grain angularity, the effective pore radii could be as low as 0.1 μm for these grain sizes, more than two orders of magnitude smaller than the resolution of the imaging technique. 13 The obtained samples measure approximately 20 mm in length and 11 mm in diameter and the image sets contain on the order of 3x108 voxels. The data are collected in the form of a stack of slices along the long axis of the sample. 11mm 20mm image slices z x y Figure 2.1. Size and orientation of core sample Using an independently determined bulk porosity value, spatial distributions of porosity can be established from the images for each voxel throughout the geologic sample. Binary porosity datasets (“binary”) are created using a threshold value above which the open pores coincide with the determined bulk porosity. For this method, porosity values in the domain assume values of 1 (open pore) or 0 (solid material). Concurrently, 16 bit porosity datasets (“porosity map”) are also created in a similar fashion. In this procedure, a porosity value ranging between 0 and 1 is assigned to each voxel. These two types of porosity datasets serve as the foundation upon which the numerical simulations are performed. 14 a) b) Figure 2.2. a) 2D image slice of 11mm diameter sandstone taken from 3D dataset obtained using MicroCT scanner. b) Section of sample showing binary spatial distribution of porosity where the white regions correspond to open pores. Ratio of white pixels to total pixels corresponds to average porosity of sample. 2.2.2 Data Cropping Given the immense amount of data acquired using the MicroCT scanner (> 108 data points), simulations can only be conducted on a subset of the porosity data due to limitations of currently available computing technology and the unavailability of a numerical model that uses parallel computing technology. Subsets of porosity data can be extracted using a cropping procedure. Cropping inevitably affects the boundary conditions for the numerical simulations, implying that mass transport across the lateral boundaries and the end of the subdomain is not possible. However, the current simulations are conceptual in nature and concentration data for comparison is not available. Cropping does therefore not affect the results obtained in this study. From the original datasets, 2 crops (subdomains) were extracted from each data set type (binary and porosity-map). 15 2.2.3 Data Binning To examine the effect of image resolution and voxel size on the simulation results, binning of the dataset is performed. Binning of the dataset is accomplished by arithmetically averaging adjacent cells in all three dimensions. If binning by a factor of two, this will halve the image resolution, double the voxel dimensions (in a 2D image), and reduce the number of voxels by a factor of 4. If starting with a binary porosity distribution as in Figure 2.2 above, the binning takes the binary data of open pores and closed media and combines the voxels together such that the data present in adjacent cells is averaged into a new larger cell. This artificially lowers the resolution of the image while subsequently increasing the number of unique porosity values to produce a more distributed porosity map as shown in Figure 2.3 below. b) a) 0 255 Figure 2.3. a) Binary porosity dataset with 144x144 voxels. Binary 8 bit histogram shows values at 0 and 255 b) Same porosity dataset after being binned 2x (to 72x72 voxels). The histogram now shows a spread in porosity values. To examine the effect of dataset resolution on the simulation results, simulations will be conducted at various resolutions based on datasets (binary and porosity map) that are binned 255 0 16 4 times by a factor of 2. Data binning is restricted to the x-y-dimension, while the original resolution is maintained along the length of the core (z-axis) (Figure 2.1) 2.3 Modeling Methods 2.3.1 Modeling Domain Tracer diffusion through the sandstone core is assessed by simulating the migration of Iodide. The transport module of the reactive transport code MIN3P is used to conduct the simulations. Using an appropriate aqueous diffusion coefficient for Iodide found in literature (1.88x10-9 m2/s, Nakashima, 2002), 3D diffusion is modeled through the subdomains of the cores. The measured porosity values are mapped onto corresponding grid nodes within the model domain. Model grid spacing depends on the resolution of the provided data. 2.3.2 Boundary and Initial Conditions Model boundaries are set such that the transport will be fully diffusion-driven. To achieve this, the hydraulic head gradient throughout the subdomain is set to zero. A constant concentration is specified at one end of the subdomain. The remaining transport boundaries of the domain are treated as impermeable no-flow boundaries. This mimics the preferred conditions for geologic burial as well as the conditions present in laboratory diffusion experiments. The initial concentration of Iodide within the domain is set to 1.0 x 10-10 mol/l and an input boundary concentration of 1.0 mol/l is specified. 2.3.3 Diffusion Modeling MIN3P uses the conservation of solute mass governing equation (equation 2.1) to calculate mass transport fluxes between the cells within the model domain. 17 QCDCV t C (2.1) Since we are only concerned with mass transport through diffusion, the velocity of the solution in our system (V) is zero. The source (Q) is kept at a constant concentration as described above. Using the experimentally acquired porosity values, along with literature values for aqueous diffusion coefficients (D0), effective diffusion coefficients (Deff) are calculated for each cell of the simulation using the Millington-Quirk approximation (Millington and Quirk, 1960) (equation 2.2): 3 4 0DDeff (2.2) Average effective diffusion coefficients are calculated across adjacent cell boundaries using either arithmetic or harmonic averaging formulations (as described in section 2.3.3). 2.3.4 Inter-Facial Fluxes - Averaging Formulations In calculating the diffusive mass fluxes between adjacent cells i and j within the simulation domain, the effective diffusion coefficients must be evaluated at the cell interfaces and may be calculated based on representative average porosities defined at the cell interface (Figure 2.4). 3 4 0, ijeffij DD (2.3) 18 φi = 1 φj = 0 di dj d φj = 0 φi = 1 iiarith dd 1 i i harm d d )( 2 1 ji 2 1 ji d 11 2 0 a) b) Figure 2.4. a) Arithmetic averaging will overestimate fluxes into areas of low porosity, b) Harmonic averaging between the cells ensures that large fluxes are not present between cells of high and low porosity Various formulations are possible to determine the average porosity. For example, the interfacial porosity can be obtained based on the spatially weighted arithmetic mean of the porosities of the two adjacent cells (equation 2.4). ij jjiiarith ij d dd (2.4) This method is commonly used in evaluating interfacial diffusion coefficients; however, as presented in Figure 2.4 for the extreme case of a fully open (water-filled) cell and a fully closed (solid-filled) cell, the average inter-cell porosity would be 50%. Calculating the effective diffusion coefficient based on equation 2.3 would yield an effective diffusion coefficient that is on the order of D0 (equation 2.5). 19 0 3 4 0 3 4 0, 39.0)5.0( DDDD ijeffij (2.5) Calculating the diffusion coefficient in such a way may overestimate diffusive fluxes into low porosity regions. Instead, the average diffusion coefficient can be calculated based on the harmonic mean (equation 2.4) j j i i ijharm ij dd d (2.6) In this case, the low porosity value dominates the harmonic average and fluxes into low porosity regions are greatly reduced. While MIN3P generally uses the arithmetic averaging formulation like many other solute transport models, the harmonic formulation has also been implemented in the model as part of the present work (Section 6.3.2). The effect of using different averaging schemes on 3D diffusion in heterogeneous media is examined through a series of separate simulations using either the arithmetic or the harmonic averaging method. Within the model, this averaging is implemented on a cell by cell basis using the diffusion coefficients at each cell as determined through the cells` porosity value and the Millington-Quirk approximation (eq. 2.2). Making use of the cell size and the distance between nodes, the model calculates the mass flux influence coefficients as equations 2.7 for the arithmetic scheme and equation 2.8 for the harmonic scheme. 20 ji j eff ji eff i ji ij xx xDxD xx A )(arith d, ji, (2.7) ieffjjeffi ji, eff j eff iharm d, ji, ΔxDΔxD ADD (2.8) As was mentioned earlier, the arithmetic averaging scheme may overestimate mass flux into low porosity regions. For this reason, the harmonic averaging scheme was favoured and was the basis for most of the simulation comparisons. For the unbinned datasets, a comparison was run between the different averaging schemes to better quantify the logical fallacies involved. 21 Chapter 3: MIN3P Diffusion Model Verification 3.1 Introduction To ensure that the model can accurately simulate diffusion-dominated transport in one, two and three spatial dimensions, it must be verified using proven methods of analysis. These generally entail the use of mathematical analytical solutions to given physical diffusion problems. Comparing the numerically derived results to the analytical solution allows for a quantification of the goodness of the fit. Also examined is the effect of the time discretization on the numerical results. The size of the time step is directly correlated to the size of the diffusion number (Steefel and MacQuarrie, 1996). When the diffusion number exceeds a value of 1, the mass transfer by diffusion during beyond the length of a single cell is significant during a single time step. This could lead to artificial numerical diffusion as is examined in section 3.3. 3.2 Verification Comparisons between the numerical diffusion model and analytical solutions are made for 1D, 2D and 3D cases. Unfortunately, while a rigorously tested analytical diffusion model can be found for the 1D case, the availability of reliable analytical diffusion models for higher dimensions is limited. When an analytical solution is unavailable, an alternative means of model verification is to compare the results to those of another widely accepted numerical model that has been previously put through its own rigorous verification. In this case, the two numerical solutions are compared when given the exact same boundary value problem. Good agreement of the solutions provides confidence that both models are solving the problem in an accurate manner. 22 3.2.1 Diffusion in 1D Semi-Infnite Domain A 1D column simulation is compared to the Ogata-Banks analytical solution (Ogata and Banks, 1961). The Ogata-Banks solution provides an analytical solution to the 1D- advection dispersion equation including linear sorption. ] 44 [),( 0 t R D t R vx erfce t R D t R vx erfcCtxC D xv (3.1) The solution is derived for a semi-infinite domain and assumes a fixed concentration boundary at the inlet, and a zero concentration initial condition. The current study focuses on diffusive transport of a non-reactive species, which allows to ignore the velocity and retardation terms. This greatly simplifies the equation to the following form: ] 4 [),( 0 Dt xerfcCtxC (3.2) Thus the concentration of the solute at a point along the column at a certain time is dependent only upon the diffusion coefficient. The parameters used for the comparison are summarized in Table 3.1. Due to the analytical codes use of bulk diffusion coefficients (Dp) as opposed to MIN3P's use of 23 aqueous diffusion coefficients (D0), conversions were made through combining equations 1.2 and 1.4 to obtain equation 3.3. 3 1 0 0 3 4 )( DDD p (3.3) These conversions were performed to ensure that the analytical and numerical models use the same effective diffusion coefficients. Output times (days) column length (m) spatial increment ∆x (m) Maximum timestep ∆t (s) Porosity (%) Diffusion coefficient (m2/s) Diffusion number 10,20,30,40,50 0.050 0.001 8640 12.5 D0 = 3.175d-11 0.274 Table 3.1. Simulation parameters used for 1D diffusion verification 24 0.0 0.2 0.4 0.6 0.8 1.0 0 0.01 0.02 0.03 0.04 0.05 Distance along Column (m) R el at iv e C on ce nt ra tio n 10 Days 20 Days 30 Days 40 Days 50 Days Figure 3.1. Diffusion at multiple output times along 1D column. MIN3P simulation (points) and Ogata- Banks analytical solution (solid lines) From Figure 3.1, we can see that the fit is very good for multiple output times between the analytical and numerical models for this simple 1D column simulation. 3.2.2 Diffusion in 2D Semi-Infinite Domain with Patch Source Initially, MIN3P results were compared to the results of PATCH3D, a 3D analytical solution (Sudicky et al., 1988). However, discrepancies between the solution and MIN3P could not be resolved and after a thorough evaluation, it was concluded that the analytical code may not yield reliable results for the parameter set chosen. As a result, the evaluation was extended to a three-way comparison and the well known groundwater flow code MODFLOW (MacDonald and Harbaugh, 1983) was used to solve the diffusion problem. 25 This approach was possible because the diffusion and groundwater flow equations are both parabolic partial differential equations, and the flow code can be used as an analogue model for the diffusion problem. 2 2 x hK t hSs Flow equation (3.4) 2 2 x cD t c p Diffusion equation (3.5) The comparison was performed by setting the hydraulic conductivity (K) equal to the bulk diffusion coefficient (Dp), and storage (Ss) is set to 1.0. Using these parameters, the calculated hydraulic head values correspond to the concentration distribution. Diffusion from a patch source into a 2D domain is simulated using MIN3P, Modflow, and Patch3D for several different diffusion coefficients. Output time (days) Column length (m) Spatial increment ∆x (m) Maximum timestep ∆t (s) Column width (m) Source Width Y (m) Diffusion coefficient (m2/s) Diffusion number D0 = 2.315d-10 0.32 D0 = 1.157d-10 0.16 D0 = 3.175d-11 0.044 50 0.10 0.0025 8640 0.40 0.03 D0 = 1.157d-11 0.016 Table 3.2. Simulation parameters used for 2D diffusion verification. 26 The model consists of a 2D strip source at x=0 with a width of 0.03 m in y-direction and is infinite in the z-direction. Diffusion occurs outwards into the x,y plane. Diffusion direction Figure 3.2. a) 2D patch source model setup. Patch source is located at x=0 and is infinite in z with a finite width in the y-direction of 7.5% the extent of the y domain. b) Diffusion contours in the x-y plane at T = 50 days x z a) y b) 27 0.0 0.2 0.4 0.6 0.8 1.0 0.0E+00 1.0E-02 2.0E-02 3.0E-02 4.0E-02 Distance along Column Centreline (m) C on ce nt ra tio n (m ol /l) D = 2.315E-10m/s D = 1.157E-10m/s D = 3.175E-11m/s D = 1.157E-11m/s Figure 3.3. Diffusion after 50 Days with multiple diffusion coefficients from a 2D patch source. MIN3P numerical simulations (points), Modflow numerical simulations (dashed lines) and Patch3D analytical solution (solid lines) For the 2D case, the three models are in good agreement when the diffusion coefficient is less than 1 x 10-10 m2/s. For larger diffusion coefficients, the analytical model Patch3D begins to provide lower solute concentrations in relation to the numerical models Modflow and MIN3P. Although at first glance this discrepancy looks to be attributable to boundary effects within the numerical models (concentration build-up at boundary), further simulations with lengthened and widened domains have discounted this as a source of error (see Appendix A). The good agreement between results obtained using the rigorously tested Modflow numerical model and MIN3P builds confidence in the numerical solutions. 28 3.2.3 Diffusion in 3D Infinite Domain with Interior Brick Source For testing the numerical code in 3D, the 3D analytical diffusion model Brick was used (Neville, 2006). For this model, a brick of finite mass is placed within the center of a system and outward diffusion in 3 dimensions is simulated. Simulations were compared for two different diffusion coefficients. Output time (days) Column dimensions x,y,z (m) Spatial increment ∆x (m) Maximum Timestep ∆t (s) Source dimensions x,y,z (m) Diffusion coefficient (m2/s) Diffusion number D0 = 1.157d-10 1.0 50 0.15, 0.15, 0.15. 0.001 8640 0.009 , 0.003 , 0.003 D0 = 3.175d-11 0.27 Table 3.3. Simulation parameters used for 3D diffusion verification. Diffusion direction x z y Figure 3.4. 3D brick source model setup. The source is located in the middle of the model domain and is in the shape of a brick with dimensions (x,y,z) of 0.009,0.003,0.003 metres. 29 Figure 3.5. The brick provides a source of finite solute mass from which diffusion occurs in all three dimensions in an expanding oblate spheroid, concentration contours shown at T = 50 days 30 0.0E+00 1.0E-03 2.0E-03 3.0E-03 4.0E-03 0.000 0.075 0.150 Length along x through brick (m) C on ce nt ra tio n (m ol /l) D = 3.175d-11 D = 1.157d-10 Figure 3.6. 3D diffusion outward from a brick source using multiple diffusion coefficients. MIN3P numerical simulations (symbols) and brick analytical solution (solid lines), T = 50 days For the 3D comparison, the diffusion of mass outwards from the finite brick source showed very similar results for both the analytical and numerical models. 3.3 Diffusion Number Analyses The diffusion number is similar in scope to the Courant number used to characterize the time scale of advective transport in relation to the grid discretization (Saaltink et al., 2001). While the Courant number is calculated using the advective linear groundwater velocity (v) (equation 3.6), the diffusion number is calculated using the systems’ pore water diffusion coefficient (Dp) (Steefel and MacQuarrie, 1996). 31 x tvCr (3.6) 2 2 x tD N pd (3.7) In both equations, and t x relate to the temporal and spatial discretization parameters respectively. For the fully implicit time integration used in MIN3P, no stability criteria exist, and simulations can be conducted with arbitrary diffusion numbers. However, large time steps may lead to artificial numerical diffusion. Conducting simulations with different diffusion numbers is an effective way to investigate the accuracy of the solution as a function of time step size. 3.3.1 1D Time Step Analysis Using the same initial parameters as in the 1D simulations presented above (Table 3.1), the model is examined using differing maximum time increments within the numerical solution. The time increments used along with their subsequent diffusion numbers are shown in Table 3.4. Maximum timestep ∆t (Days) 0.0005 0.001 0.005 0.01 0.05 0.1 0.5 1 5 Diffusion number 0.0014 0.0027 0.014 0.027 0.14 0.27 1.4 2.7 14 Table 3.4. Simulation parameters used for 1D diffusion verification Figure 3.7. depicts the percentage error between the simulation with the smallest maximum time increment and the subsequent simulations with larger time increments. The results 32 demonstrate that relatively large percentage errors occur near the leading edge of the diffusion front for simulations with large diffusion numbers. However, given the magnitude of the concentrations involved at the diffusion front, these discrepancies can be considered to have very little effect on the overall transport of the solute and the simulation results. Comparing Figures 3.7 and 3.8, this correlation between concentration and percentage error can be seen as the percentage errors only increase significantly in the second half of the column, where the solute concentration is lowest. -50 0 50 100 150 200 250 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 Distance along Column (m) % E rr or a s co m pa re d to ti m es te p 0. 00 05 0.001 Days 0.005 Days 0.01 Days 0.05 Days 0.1 Days 0.5 Days 1 Day 5 Days Figure 3.7. Percent error after 50 days along 1D column. MIN3P simulations using differing maximum time steps. 33 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 Distance along Column (m) N or m al iz ed c on ce nt ra tio n 0.001 Days 5 Days Figure 3.8. Concentration profile after 50 days along 1D column. MIN3P simulations using differing maximum time steps. 3.3.2 3D Time Step Analysis Results for the simple 1D homogeneous case are encouraging, but cannot be extrapolated to 3D heterogeneous media. Using a 3D dataset with heterogeneous porosity distribution, an additional analysis was completed to determine the effect of the diffusion number on 3D concentration distributions. This test is directly related to the main focus of this study due to the fact that the material examined in this time-step analysis is the same dataset that is used in the remainder of this study. A smaller dataset that had been derived from the original (section 4.2.1 – cropping and binning) was chosen so that simulations with 34 small minimum time increments could be completed without need for long computation times. Binned amount # of cells (x by z) cell size (m) bulk porosity (%) bulk effective D (m2/s) Maximum time increment (Hours) Diffusion number 0.0001 0.159 0.0005 0.79 0.001 1.59 4x 36 by 24 7.20E-05 22.3974 2.57E-10 0.01 15.9 Table 3.5. Simulation parameters used for 1D diffusion verification The porosity and subsequent effective diffusion coefficient are determined through the mean value of the binned dataset. Due to its origin as a binary dataset, the full range of porosities between 0 and 1 are actually present within the column. Subsequently, the plotted errors will be exacerbated in areas of connected high porosity where the effective diffusion coefficient can be as high as 1.88x10-9 m2/s As with earlier time step analyses, a diffusion profile is simulated using MIN3P with several differing maximum time increments. These are examined by comparing the results at each time increments to those of the smallest time increment (Figure 3.9). A contour plot comparing 2D concentration profiles after 0.1 hours is also examined (Figure 3.10). For the 3D simulations in subsequent sections, a maximum time increment of 0.1 hours was used to ensure reasonable computation times. Given the variable nature of MIN3P’s time integration scheme, this maximum value was reached only after the main diffusion event examined in the results had occurred. 35 -3% -2% -1% 0% 1% 2% 3% 4% 5% 6% 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (m) % Er ro r a s co m pa re d to ti m es te p of 0 .0 00 1h rs 0.0005Hrs 0.001Hrs 0.01Hrs Figure 3.9. Percent error after 0.1 Hours along a 3D column. MIN3P simulations using differing maximum time step. Figure 3.10. Contour plot of concentrations after 0.1 hours using time steps of 0.0001 (black lines) and 0.01 hours (white lines). Cross section is through the middle of the column with source located at z = 0 m. 36 The effect of the increasing diffusion number can be seen in Figure 3.9 as the discrepancies grow larger with each increasing time increments. As was determined with the 1D case, the largest relative discrepancies occur near the end of the column where concentrations are at their lowest and can be considered to have very little effect on the overall transport of the solute and the simulation results. Comparing the black and the white contours in Figure 3.10, the most visible errors occur where the area through which the solute is diffusing is at its largest. This corresponds to large areas of higher porosity. These significant errors can occur within the body of the diffusion in the first half of the column. Through interpolation of the averaged results above, it is recommended that a diffusion number of 10 not be exceeded if the percent error is to be kept at less than 1%. 37 Chapter 4: Sub-mm Scale Diffusion Modeling in a Sandstone Core 4.1 Introduction Due to its relative large grain size, a sandstone core was chosen as the medium of interest for this study. The 11mm x 20mm core was examined at UNB through non- destructive microCT imaging to produce high-resolution images of the internal porosity structure of the sample. The raw dataset produced by the microCT imager was processed using knowledge of the medium’s bulk porosity. Combining the raw image data with the bulk porosity produced two types of porosity datasets which could be used as the basis for the simulations presented in the following. After receiving the datasets from UNB, they were processed further through cropping and binning into a series of 10 separate datasets for each type (binary and porosity map). Each of these sub datasets was subjected to a number of simulations. A summary of the completed simulations is presented in Table 4.1. 38 Type Crop Location Binned amount # of cells (x by z) Bit Depth bulk porosity (%) Binary Crop1 1x 144 by 96 1 22.5266 2x 72 by 48 3 22.4814 4x 36 by 24 6 22.3974 8x 18 by 12 9 22.2950 16x 9 by 6 12 22.2300 Crop2 1x 144 by 96 1 10.1696 2x 72 by 48 3 10.1434 4x 36 by 24 6 10.0862 8x 18 by 12 9 9.9854 16x 9 by 6 12 9.8729 PorosityMap Crop1 1x 144 by 96 16 16.3521 2x 72 by 48 19 16.3513 4x 36 by 24 22 16.3506 8x 18 by 12 25 16.3500 16x 9 by 6 28 16.3495 Crop2 1x 144 by 96 16 11.3812 2x 72 by 48 19 11.3805 4x 36 by 24 22 11.3798 8x 18 by 12 25 11.3793 16x 9 by 6 28 11.3787 Table 4.1. Summary of porosity data sets created for simulations The bit depth refers to the number of unique possible porosity values that the model contains. The actual maximum number of unique values can be determined by calculating 2bit depth. For example, a bit depth of 3 has 23 = 8 unique values for porosity. 4.2 Dataset Description Processing of the raw microCT data at UNB ended up in two separate datasets. While the resolution of the microCT was not sufficient to resolve individual grains and pores, pore clusters are resolvable and produce characteristic response values in the raw dataset. This dataset was processed using a threshold value for porosity above which the voxel was defined as an open space (i.e. void), and the remaining voxels were defined as solids. The 39 threshold was chosen such that the relative volume occupied by voids corresponds to the bulk porosity value. This approach is common in image processing and has previously been used by various authors in porous media simulations (Nakashima, 2002, Al-Raoush and Willson 2005). For the 16bit porosity map dataset, the threshold was lowered and the voxels above the threshold were given a value of porosity based on their raw image value (grayscale value). These were again scaled to match the bulk porosity of the sample. This gives a dataset with raw data values representative of the porosity distribution that range between 0 and 65535. These values can subsequently be normalized to between 0 and 1 (65535 values between 0 and 1). From these datasets, 2 crops were extracted and subsequently binned. 4.2.1 Cropping and Binning The original dataset consists of 1x108 voxels (379x379x799) and is too large to be directly used in diffusion simulations. Two separate sub-datasets of 144x144x96 voxels (2x106) each are extracted from the full dataset to facilitate numerical analysis (Figure 2.2). 40 Figure 4.1. Two separate 2x106 voxel datasets cropped from the original 1x108 voxel MicroCT porosity dataset. The top and bottom binary datasets have overall average porosities of 10% and 22% respectively. The crops (subdomains) are located 20 pixels away from the edges of the full dataset that they are adjacent to. This choice was made to avoid possible boundary effects of the imaging techniques, which may affect the quality of the porosity data. Crop 1 is located in the lower right hand corner, while Crop 2 is in the upper left. While both the full binary and 16 bit datasets have comparable bulk porosities of approximately 14.5%, the binary dataset shows much more variability between the sub-domains. The average bulk porosity of the slices of Crop 1 is 22%, while the average bulk porosity of the slices of Crop 2 is only 10%. A similar comparison of the crops of the 16 bit dataset (porosity map) only shows average bulk porosity values of 16% and 11 %, respectively. The smaller variance of the 16 bit dataset is due in part to resolution limitations in the imaging technologies which decreases the effect of any averaging scheme emplaced on the data. The subdomains were purposefully chosen to capture a relatively large difference in bulk porosity. Due to the appearance of edge effects encountered near the first and last slices of the sample caused by imaging 41 inaccuracies, the 96 slices in the long (z) direction were taken starting from the 200th slice from the bottom. Using the process explained earlier, each dataset is binned repeatedly until a resolution is reached that is too coarse to be modeled effectively. For the datasets in question, this occurs after 4 binning iterations when the original dataset has been binned by a factor of 16. This degree of binning leads to coarse datasets containing 9x9x6 voxels. As will be examined later, although some parameters remain relatively unchanged with the degree of binning, the loss of resolution can have a significant impact on investigations into the transport processes. Irrespective of the data format and amount of binning, the modeling method remains the same (section 2.3). The cropping and binning of the dataset in the X and Y dimensions was completed using built-in tools in the image manipulation program imageJ (http://rsbweb.nih.gov/ij/). In the Z dimension, a program developed for this project and outlined in Appendix B.2 was used to bin the separate slices. 4.2.2 Binary Datasets The binary dataset loosely represents the pore structure of the sandstone media. It is not a perfect representation due to the resolution limitations of the microCT device. For a true pore-scale system to be represented, the resolution of the imaging device has to be high enough such that the individual pores can be resolved. For diagnostic purposes, the binary dataset does represent a system of open and closed pores, but it may not be representative of our sample. 42 0 20 40 60 80 100 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column, z (m) Po ro si ty (% ) 0 20 40 60 80 100 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column, z (m) Po ro si ty (% ) a) b) Figure 4.2. Image of porosity of first slice and averaged porosities per slice versus distance along column. a) Crop 1 of the binary dataset. Average porosity of 22% b) Crop 2 of the binary dataset. Average porosity of 10%. Once cropped, the two datasets are each binned 4 times by a factor of 2. Each binning decreases the resolution in each of the coordinate directions (x,y,z) by 2 and consequently increases the voxel volumes by a factor of 8. Thus the total number of voxels lost due to the averaging during each binning iteration is 7/8th of the dataset prior to binning. Binned amount # of cells (x by z) cell size (m) bulk porosity Crop 1 (%) bulk porosity Crop 2 (%) 1x 144 by 96 1.80E-05 22.53 10.1696 2x 72 by 48 3.60E-05 22.48 10.1434 4x 36 by 24 7.20E-05 22.40 10.0862 8x 18 by 12 1.44E-04 22.30 9.9854 16x 9 by 6 2.88E-04 22.23 9.8729 Table 4.2. Binning of binary dataset crop 1 43 4.2.3 16 Bit (Porosity Map) Datasets The 16 bit dataset attempts to represent the porosity distribution of the material by assigning a range of porosity values to the images, again based on the bulk porosity. It assumes that the imaging technique is not precise enough to resolve individual grains and pores and thus weighs the actual porosity of the individual voxels according to the raw value as determined using the imaging techniques. Cropping this dataset into Crop 1 and Crop 2 as with the binary dataset, the effect of the additional 65535 bits of data become apparent as the variance of the dataset is greatly reduced. The two crops have respective average porosities of 16.5% and 13.5% while the original had a bulk porosity of 15% 0 20 40 60 80 100 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column, z (m) Po ro si ty (% ) 0 20 40 60 80 100 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column, z (m) Po ro si ty (% ) a) b) Figure 4.3. Image of porosity of first slice and averaged porosities per slice versus distance along column. a) Crop 1 of the 16 bit dataset. Average porosity of 16.5% b) Crop 2 of the 16 bit dataset. Average porosity of 13.5%. 44 Like the binary dataset, once cropped, the datasets are each binned 4 times by a factor of 2. Binned amount # of cells (x by z) cell size (m) bulk porosity Crop 1 (%) bulk porosity Crop 2 (%) 1x 144 by 96 1.80E-05 16.3521 11.3812 2x 72 by 48 3.60E-05 16.3513 11.3805 4x 36 by 24 7.20E-05 16.3506 11.3798 8x 18 by 12 1.44E-04 16.3500 11.3793 16x 9 by 6 2.88E-04 16.3495 11.3787 Table 4.3. Binning of 16 bit dataset crop 1 4.3 Diffusion modeling Once cropped and binned, the datasets are formatted such that they can be read into the reactive transport code MIN3P. This formatting is done using a program developed for this project and outlined in Appendix A.2. MIN3P is used for the simulations because it allows the inclusion of multiple components and reactions, which may be used in future studies. However, the present work only focuses on non-reactive diffusive solute transport. Using Iodide as a conservative tracer, simulations are created to conduct numerical diffusion experiments. Two different set-ups are investigated, including a single reservoir and a double reservoir (through-diffusion) design. The single reservoir system simulates a closed system from which Iodide can only enter from the input end and will eventually build up in the system until it is at a constant concentration throughout the medium. Simulations of through-diffusion provide a numerical analogue to more complex double reservoir experiments. With a second reservoir present at the end of the domain, the system is allowed to reach a steady state much sooner where the inputted Iodide mass is 45 balanced by an equal amount of mass exiting into the low concentration reservoir. This type of system allows for better comparisons of mass flux between the differing model media. Starting with the medium’s pore structure saturated with water, the reservoirs are kept at a constant concentration of Iodide along one face of the sample. The Iodide is allowed to diffuse into the medium based on the laws that govern its movement constrained by the inputted parameters and the model criteria. 144x144 cells Constant Iodide concentration of 1E-10 mol/l on top face No mass flux through vertical sides 96 cells z Constant Iodide concentration of 1.0 mol/l on bottom face x y Figure 4.4. Model domain setup for numerical through diffusion experiment with Iodide reservoir input at z = 0. 4.3.1 Input Parameters Along with the porosity maps generated from the imaging techniques, the diffusion model is influenced by several other input parameters. The aqueous diffusion coefficient of Iodide is the diffusion coefficient of Iodide in pure water (no solids. porosity of 1.) It has been 46 calculated and tabulated in prior laboratory experiments (Nakashima, 2002). Effective diffusion coefficients are calculated based on the Millington relationship as outlined in Section 2.3.3. Advective transport was excluded by setting the head gradient across the domain to zero and specifying low hydraulic conductivity for the media (10-10 m/s) The boundary conditions of the domain consist of zero flux boundaries for each of the vertical faces (sides of the domain). For the through-diffusion experiments, constant concentration boundaries are located at the top and bottom face, the bottom face at z = 0 being coincident with a constant Iodide concentration of 1 mol/l. Paralleling a flow-through diffusion cell experiment, the top face is kept at a negligible concentration such that any Iodide reaching the end of the domain will exit the system. For the closed system formulation, a constant concentration boundary is only specified at the inlet while the end of the domain is assumed impermeable (zero concentration gradient). The simulation is run until a final solution time of 1.0 hour with intermediate output times of 0.005, 0.01, 0.02, 0.04, 0.1, and 0.5 hours. This allows for examinations of contaminant movement at early diffusion times while also capturing the system at quasi steady state conditions. Longer solution times would better show the steady state diffusion transport, but were not attempted due to constrains on model runtimes. Table 4.3 provides a summary of the model input parameters. Simulation output times (hours) Simulation final output time (hours) Diffusion coefficient D0 (m2/s) Head gradient across domain (m/m) Input concentration z = 0 (mol/l) Output concentration z = max (mol/l) 0.005, 0.01, 0.02, 0.04, 0.1, 0.5, 1.0 1.0 1.88E-09 Zero 1.0 1.0E-10 Table 4.4. Major input parameters for sandstone core diffusion simulations 47 The maximum simulation times and output times were chosen after conducting numerous test simulations to provide information on the rate of diffusion and the duration required to achieve quasi-steady state conditions (for through-diffusion experiments). The output times were chosen such that the evolution of the diffusing Iodide could be examined as it entered the media. For the closed system, the final run time was chosen such that the amount of Iodide reaching the end of the domain would be minimal. In the open system, the chosen final run time allowed for a number of the faster simulations to approach steady state. Longer runtimes were also avoided to minimize any constraints due to lengthy computational runtimes. 4.4 Results For the single sandstone core of 11 mm by 20 mm, a large number of simulations were completed. With multiple crops of both the binary and 16 bit porosity datasets, changes in the degree of binning as well as variations of the model workings allow for the simulations to be compared. 4.4.1 Simulations of 3D Diffusion for Binary Dataset Using the provided binary datasets as porosity maps, simulations of diffusion through the datasets were completed for diffusion timelines of up to 1 hour. Given the binary nature of the dataset’s porosity, diffusion through the sample was not expected to occur along a planar diffusion front as would be the case in a homogeneous system, but instead be determined by the pore structure itself. To examine the movement of the Iodide as it diffuses through the sample, a slice through the sample is examined through time and a number of slices through the sample are 48 examined at the end of the simulation. This is done for both the binary sample dataset as well as the sample dataset that has been binned by a factor of 4. Figure 4.5 shows the changes in the concentration of Iodide in one slice 3.3x10-4 m from the Iodide source as the Iodide diffuses into the sample through time. Figure 4.5. Concentration map at xy slice 19 from the source (z=3.3x10-4m) after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the binary dataset. These images show the heterogeneity of the diffusional flux through the sample’s pore structure. The warmer colours of higher concentration are not always found in most open areas. As such, the concentration does not seem to depend on the porosity of an area alone, but also the connectivity of that area to pores above, below, and behind it. This can be further examined through looking at the slices along the transport path through time. Figure 4.6 shows the changes in the concentration of Iodide along the centreline of the sample in the y direction as the Iodide diffuses into the sample through time. 49 Figure 4.6. Concentration map at xz slice along the centerline of y after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the binary dataset. These images along the transport path greatly depict the three dimensionality of diffusive transport, even for a fairly homogeneous sandstone. By looking at only one plane of cells, we can see that the Iodide, entering from the bottom, does not enter homogeneously all along the plane and does not stay in the chosen plane as it moves upwards through the sample. This is due to the movement perpendicular to the plane as the diffusing molecules follow the tortuous paths of the areas of connected porosity. Looking at multiple slices through the sample after diffusion is allowed to occur for 1.0 hours (Figure 4.7), we again can see that the areas of warmest colour, and thus highest concentration, vary from slice to slice. 50 Figure 4.7. Concentration map after 1hr through xy slices a) 13 b) 26 c) 39 d) 52 e) 65 f) 78 and g) 91 out of 96 for the binary dataset. As the concentration decreases with distance, the areas of high concentration migrate around the slices from the lower left in the slice near the input to the upper left in slices near the middle of the sample to the centre right in slices farthest away from the source, indicating the presence of preferential diffusion pathways. If these preferential pathways did not exist, the Iodide would take the shortest route away from the source and would thus not move in such a way around the sample. 51 Figure 4.8. Concentration map at xy slice 19 from the source (z=3.3x10-4m) after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the dataset that has been binned by 4x. When the dataset is binned by a factor of 4 as in these images, the percentage of the slice that is available for Iodide penetration is greatly increased. Even so, the Iodide diffusion front still does not move through the slice all at once. There are several areas where faster moving fingers of Iodide penetrate and move through the slice before laterally adjacent pores are filled. This displays the continued reliance of the diffusion path on the connectivity of the pores. Again, this can be further examined through looking at the slices along the transport path through time (Figure 4.8). 52 Figure 4.9. Concentration map at xz slice along the centerline of y after a) 0.005hr b) 0.01hr c) 0.02hr d) 0.04hr e) 0.1hr and f) 0.5hr for the dataset that has been binned 4x. In these xz slices along the y centreline, the Iodide is present throughout a larger area of the cross section when compared to the similar un-binned case (Figure 4.6), but there is still the presence of fingering at the diffusion front where the Iodide preferentially follows the paths of least resistance. This connectivity of the pores in 3-dimensions can be overlooked in many groundwater problems where the REV is chosen such that it does not encompass these heterogeneities. But for small scale transport such as this study, these preferential pathways are the main conduits for diffusion through the sample and it is even probable that these paths at these small scales will make the diffusion anisotropic. With multiple crops of both the binary and 16 bit porosity datasets, changes in the degree of binning as well as variations of the model workings all have an effect on the diffusion of Iodide through the sample. The subsequent sections examine and compare the various simulation parameters and their effects. 53 4.4.2 Crop 1 Crop 1 is located in the lower right quadrant of the original dataset. It is an area of higher than average porosity with denser clusters of open pore spaces. 4.4.2.1 Binary Dataset The Crop 1 binary dataset has the highest average porosity of all of the datasets at 22%. All things being equal, it should potentially show the highest rate of diffusion. Binned amount # of cells (x by z) cell size (m) bulk porosity (%) Runtime Cumulative Mass In (mol/1.0hr) 1x 144 by 96 1.80E-05 22.5266 39 hrs 37 min 12 sec 1.87E-06 2x 72 by 48 3.60E-05 22.4814 1 hrs 0 min 44 sec 2.19E-06 4x 36 by 24 7.20E-05 22.3974 2 min 46 sec 2.68E-06 8x 18 by 12 1.44E-04 22.2950 11.39 sec 3.23E-06 16x 9 by 6 2.88E-04 22.2300 0.83 sec 3.73E-06 Table 4.5. Binned parameters for binary dataset crop 1 Figure 4.10 shows the changes in the concentration distribution of Iodide due to the binning of the dataset. Figure 4.11 shows how the concentration along the long axis of the domain is affected. As the dataset becomes more binned, the number of unique porosity values increases while the total number of voxels decreases. With the lowering spatial resolution, the amount of data available to the model decreases correspondingly (number of voxels) resulting in a less complex model and shorter runtimes. In order to maintain accuracy, there is more information held within each voxel due to their combination while binning (increase bit depth). 54 Figure 4.10. Concentration map at an xy slice 3.3x10-4m from the source after 0.01hrs for the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x. As the number of unique porosity values in the input database increases, the volume of zero-porosity decreases and the connectivity of the open pores increases. This results in the spreading of the solute through larger areas of the sample as seen in Figure 4.10. This also results in faster diffusion as is visible in the increases in the mass transfer rates in Table 4.5 above as well as in the flattening of the curves in Figure 4.11 below. Due to the solute transport in the column being confined to diffusion, these increases in mass transfer rates correspond to an increase in the bulk effective diffusion coefficient. Since the bulk porosity remains relatively constant throughout the simulations, the increase in the mass transfer rates must be dependent on an increase of connectivity between the open pores within the column. 55 00.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour Figure 4.11. Concentration profiles along the long axis of the domain for each output time for the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x 4.4.2.2 16 Bit Dataset The Crop 1 16 bit dataset, while in the same location as the binary dataset, only has a bulk porosity of 16%. This is higher than the bulk porosity of the full dataset`s 14%, but still shows nowhere near the variance presented in the binary dataset. 56 Binned amount # of cells (x by z) cell size (m) bulk porosity (%) Runtime Cumulative Mass In (mol/1.0hr) 1x 144 by 96 1.80E-05 16.3521 20 hrs 38 min 55 sec 2.68E-06 2x 72 by 48 3.60E-05 16.3513 45 min 9 sec 2.66E-06 4x 36 by 24 7.20E-05 16.3506 2 min 23 sec 2.67E-06 8x 18 by 12 1.44E-04 16.3500 10.11 sec 2.71E-06 16x 9 by 6 2.88E-04 16.3495 0.88 sec 2.75E-06 Table 4.6. Binned parameters for 16 bit dataset crop 1 Since the action of binning is to arithmetically average adjoining cells, the differences seen in the binary dataset are not visible here. Due to its higher starting bit depth, the 16 bit dataset is not as affected by the increase in internal voxel information as was the case in the binary dataset. This leads to less overall changes within each voxel due to their combination as the total number of voxels decreases. This is also displayed in the relatively constant values of the cumulative mass in. The porosity map information as determined by the spatial position of the voxels is not influenced by the binning procedure and the solute diffusion into the system remains relatively unaffected. 57 Figure 4.12. Concentration map at an xy slice 3.3x10-4m from the source after 0.01hrs for the porosity dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x. As with the binary dataset, the binning affects the resolution of the datasets and leads to a smaller amount of voxels. Contrastingly, the pores are spread relatively uniformly through the dataset, resulting in a limited volume of zero-porosity through which the solute does not travel. This results in a fairly constant mass flux into and through the sample as is visible in Figures 4.12 and 4.13. 58 00.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour Figure 4.13. Concentration profiles along the long axis of the domain for each output time for the porosity map dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x There is very little variation between the diffusion profiles through the datasets binned at different amounts. This agrees with the relatively constant cumulative mass fluxes from Table 4.6 as well as the connectivity of the pores in the system as examined through the ratios in Table 4.7 below. 59 4.4.2.3 Binary Versus 16 Bit Datasets The two datasets create quite different models for diffusion. With its pore structure made up of either open pore space or closed solid media, the unbinned binary dataset forces the Iodide diffusion to travel only through areas of connected pores. Transport through the solid media is minimized and hence somewhat convoluted pathways may emerge. On the other hand, the amount of purely solid media present in the 16 bit dataset is of a much smaller proportion. Much of the solid matrix that is found in the binary dataset is instead already averaged alongside the open pore structure. This pore structure, while having a lower porosity than in the binary case, is more even throughout the media. Thus the transport path for Iodide diffusion is less through specific pore pathways and occurs more as a unified diffusion front. This difference disappears as the datasets are binned and the diffusion front model becomes prevalent throughout. This can be seen in the ratios calculated in Table 4.7 below. The ratio of Porosity to Cumulative Mass In gives insight into how easily the solute can enter and travel along the open pores of the system. The higher ratio signifies that the pores are more connected and that it is relatively easy for the solute to enter and move about the system. Binary 16 Bit Binned amount bulk porosity (%) Cumulative Mass In (mol/1.0hr) Ratio (Porosity / mass in) bulk porosity (%) Cumulative Mass In (mol/1.0hr) Ratio (Porosity / mass in) 1x 22.5266 1.87E-06 8.3E-6 16.3521 2.68E-06 16.4E-6 2x 22.4814 2.19E-06 9.7E-6 16.3513 2.66E-06 16.3E-6 4x 22.3974 2.68E-06 12.0E-6 16.3506 2.67E-06 16.3E-6 8x 22.2950 3.23E-06 14.5E-6 16.3500 2.71E-06 16.6E-6 16x 22.2300 3.73E-06 16.8E-6 16.3495 2.75E-06 16.8E-6 Table 4.7. Comparison of ratios of porosity / mass in From Table 4.7, the unbinned binary dataset has the lowest ratio and thus the lowest amount of interconnectivity between its pores. This occurs due to the large volume of zero 60 porosity material within the dataset. As the binary dataset is binned, the volume of zero porosity material decreases through the averaging effect of the binning and the ratio as well as the connectivity increase. Within the 16 bit porosity map, there is very little variation between the ratios of the binned datasets. This occurs because the distribution of porosity within the system is initially well spread and thus subsequent binning has little effect on the connectivity of the pores. 4.4.2.4 Arithmetic Versus Harmonic Averaging Similar spreading of diffusion into areas of lower porosity can be seen when comparing the arithmetic and harmonic averaging schemes. For the binary case, the difference between the averaging schemes is very noticeable as the connectivity of the separate pores is artificially increased by diffusion through areas of zero-porosity. Dataset Type Dataset Type bulk porosity (%) Runtime Cumulative Mass In (mol/1.0hr) Harmonic 22.5266 39 hrs 37 min 12 sec 1.87E-06 Binary Arithmetic 22.5266 28 hrs 37 min 25 sec 2.62E-06 Harmonic 16.3521 20 hrs 38 min 55 sec 2.68E-06 16 bit Arithmetic 16.3521 25hrs 39 min 51 sec 2.76E-06 Table 4.8. Parameters of arithmetic vs harmonic Due to the binary nature of the dataset, areas of open pore are frequently located adjacent to areas of zero porosity matrix. As explained earlier, arithmetic averaging leads to artificial diffusion into the areas of zero or low porosity. This can be seen in the comparisons in Figure 4.14 as the slice from the model using arithmetic averaging shows large areas with high Iodide concentration in regions that have low porosity. The Iodide diffusion path from 61 the model using harmonic averaging on the other hand more closely follows the apparent pore structure of the media. Figure 4.14. Slice of porosity 3.3x10-4 m from the source for binary data set, as well as concentration maps for models with harmonic and arithmetically averaged modeling schemes (respectively) For the 16 bit case, the difference between the arithmetic and the harmonic averaging schemes is not as prevalent as in the binary case. Given that the adjacent cells are more likely to be similar in value in this dataset than in the binary case, the choice of averaging scheme has a greatly reduced effect. This could be an artifact of the original processing step used to create the 16 bit dataset from the raw microCT data. Due to the resolution limits emplaced by the imaging device, the dataset has already undergone a certain amount of “averaging”. Figure 4.15. Slice of porosity for 16 bit data set (porosity map), as well as concentration maps for models with harmonic and arithmetically averaged modeling schemes (respectively) 62 This inferred imaging device limited averaging can be compared to the arithmetic averaging of the porosity maps that occurs through the binning process. While it’s argued above that harmonic averaging is best for contending with artificial diffusion into areas of zero or low porosity during solute transport, harmonic averaging during the binning process would quickly close off the open pores and lower the overall bulk porosity of the sample as its connectivity is reduced to zero. Arithmetic averaging, conversely, keeps the bulk porosity of the sample relatively constant throughout the binning process, but it also spreads the open pores among a greater volume of the sample thus potentially increasing the connectivity in the process. 4.4.3 Crop 2 Crop 2 is located in the upper left quadrant of the original dataset. It is an area of lower than average porosity with small areas of connected pores. 4.4.3.1 Binary Dataset The Crop 2 binary dataset has the lowest bulk porosity value of all of the datasets with a porosity of 10%. With the low bulk porosity and thus the large areas of zero porosity solid matrix, the unbinned binary dataset has a very restricted transport pathway through the 10% of the media that is open to diffusion (Figure 4.10). 63 Binned amount # of cells (x by z) cell size (m) bulk porosity (%) Runtime Cumulative Mass In (mol/1.0hr) 1x 144 by 96 1.80E-05 10.1696 65 hrs 21 min 50 sec 2.13E-07 2x 72 by 48 3.60E-05 10.1434 1 hrs 41 min 40 sec 3.20E-07 4x 36 by 24 7.20E-05 10.0862 3 min 22 sec 5.22E-07 8x 18 by 12 1.44E-04 9.9854 12.06 sec 8.17E-07 16x 9 by 6 2.88E-04 9.8729 0.98 sec 1.12E-06 Table 4.9. Parameters for binary dataset crop 2 as a function of binning degree Similar to Crop 1, as the dataset becomes more binned, the volume of voxels of zero- porosity decreases as the binning process arithmetically averages them with surrounding open pores. Due to its low starting porosity, the number of pores open to influx of solute from the source is relatively small. This can be seen in the small value for cumulative mass in. Once binned 16x, the area open to influx of solute from the source is greatly increased and the cumulative mass influx into the system increases by a factor of 5. Figure 4.16. Concentration map for crop 2 at an xy slice 3.3x10-4m from the source after 0.01hrs for crop 2 of the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x. 64 Through arithmetic averaging, the binning quickly appears to increase the area through which the Iodide can diffuse. Seen clearly in the concentration profile (Figure 4.16), each subsequent binning greatly increases the transport of the Iodide (Table 4.8). Though it may seem detrimental, arithmetic averaging is used such that the bulk porosity remains relatively constant throughout the binning. 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour 0 0.2 0.4 0.6 0.8 1 0.0E+00 4.0E-04 8.0E-04 1.2E-03 1.6E-03 Distance along column (z) N or m al iz ed C on ce nt ra tio n 0.005 Hours 0.01 Hours 0.02 Hours 0.04 Hours 0.1 Hours 0.5 Hours 1 Hour Figure 4.17. Concentration profiles for crop 2 along the long axis of the domain for each output time for the binary dataset. Binned a) 1x b) 2x c) 4x d) 8x and e) 16x. 65 4.4.3.2 16 Bit Dataset With a bulk porosity of 12%, Crop 2 again has a lower porosity than the average of the original dataset. Due to the smaller variance of the 16 bit datasets, the lower porosity of Crop 2 acts similar to its higher porosity counterpart. Like the Crop 1, subsequent binning does not have much of an effect on the porosity and thus on the Iodide transport. Binned amount # of cells (x by z) cell size (m) bulk porosity (%) Runtime Cumulative Mass In (mol/1.0hr) bin1x 144 by 96 1.80E-05 11.3812 22 hrs 3 min 44 sec 1.52E-06 bin2x 72 by 48 3.60E-05 11.3805 39 min 43 sec 1.51E-06 bin4x 36 by 24 7.20E-05 11.3798 2 min 31 sec 1.53E-06 bin8x 18 by 12 1.44E-04 11.3793 10.62 sec 1.59E-06 bin16x 9 by 6 2.88E-04 11.3787 0.89 sec 1.65E-06 Table 4.10. Binned parameter for 16 bit dataset crop 2 4.4.3.3 Binary Versus 16 Bit With the lower bulk porosity of Crop 2, quite a noticeable difference can be seen between the 10% porosity binary dataset and the 12% porosity 16 bit dataset. While the large amount of solid matrix present in the binary dataset greatly limits the transport through the media, the slightly lowered average porosity of the 16 bit dataset has little effect on lowering the rate of the Iodide transport through its domain. The contrast between the two dataset types can be attributed to the differences in the connectivity of their open pores. The consequence of this difference depends greatly on the resolution of the sample in question and the scale that is of interest. The binary dataset needs high resolution imaging to be sure that thresholding can adequately describe the pore structure of the sample. Due to the large computational constraints required by high resolution images, this limits the overall size 66 of the sample which in itself can create experimental challenges. To bypass some of these challenges, Nakashima (2003, 2008) chose to examine the diffusion of small numbers of particles within certain samples using random walk simulations. While this effectively characterizes the diffusion characteristics of the rocks, it lacks the ability to be expanded beyond the realm of small scale solute transport. The 16 bit dataset doesn’t depend as much on the resolution of the images since each voxel is already a mixture of open and closed pore space. From the perspective of larger scale averaging and the effectiveness of binning, the 16 bit dataset appears to have the advantage of robustness and consistency. But, from a small-scale connectivity perspective, a binary dataset appears to better maintain tortuous transport pathways which are lost in the averaging of the binned as well as the 16 bit datasets. Unfortunately, the simulated behaviour of the solute through the sample has yet to be experimentally examined and as such it is unclear as to how much the connectivity of the pores influences the diffusion pathways at the scale of this project. With experimentally derived solute concentrations, this could be determined. Unfortunately, these were not available at this time. 67 4.4.4 Crop 1 Versus Crop 2 The bulk porosity of Crop 1 for the binary and 16 bit datasets is 20% and 16% respectively. This differs from the values of 10% and 12% found for Crop 2. This difference is much more visible on the transport within the binary datasets. The lack of connected pathways of open pores in Crop 2 substantially slows the progress of Iodide diffusion through the media. Contrastingly, within the 16 bit datasets, both crops have a lack of areas of zero porosity and thus the transport is not limited to any specific areas and the Iodide transport is fairly even throughout the media for both crops. Crop Dataset Type bulk porosity (%) Runtime Cumulative Mass In (mol/1.0hr) Binary 22.5266 39 hrs 37 min 12 sec 1.87E-06 1 16 bit 16.3521 20 hrs 38 min 55 sec 2.68E-06 Binary 10.1696 65 hrs 21 min 50 sec 2.13E-07 2 16 bit 11.3812 22hrs 3 min 44 sec 1.52E-06 Table 4.11. Comparisons between crops for unbinned datasets 68 4.5 Discussion Using Fick's law to model diffusion, studies have been able to adequately mimic the bulk diffusive transport visible in bench and field scale experiments (Voss and Souza, 1987. Zhang et al, 1998. Harvey and Gorelick, 2000, Patriarche et al., 2004.). When looking at the pore scale, further success has been found by using random walk simulations in place of Fick's law (Nakashima, 2003). Given the current state of easily available computational capacity, expansions of pore scale simulations are hitting a ceiling. Thus, in trying to fill the gap, this study has examined bringing Fick's law to the local REV scale (porosity-map dataset). The unbinned binary dataset simulations essentially correspond to the application of Fick’s law to the pore scale. To ensure that the pore scale model adequately represents the samples in question, a pore size analysis must be conducted to make sure that the resolution of the investigating instruments is high enough that it can resolve the individual pores. If this is not the case, a binary dataset can be constructed, but it will not represent the sample correctly. With poorly resolved pore/matrix boundaries, the distinction between open and closed pores becomes unclear. The thresholding used to create a binary dataset can then create smaller or larger pores dependent on where the threshold is set. This can affect the connectivity of the sample as pathways can be either closed off or enlarged respectively. By using Fick's law, the main driving force for diffusion is the concentration gradient that is set up within the media. While the binary dataset may not be pure pore scale modeling due to limitations in the resolution of the investigative techniques, assessments can be made of the diffusive transport mechanisms of contaminants on the scale of the materials' pore structure. Given the transportation timelines and velocities involved, this is the scale of interest. 69 Unfortunately, without pertinent concentration profile data, the assessments made cannot be directly compared to real-world laboratory results and thus cannot at this time be verified for correctness. Once this type of data can be reliably obtained, the errors inherent in each of the tested modeling schemes can be better quantified and compared to actual transport characteristics. Even though this REV scale modeling is less computationally intensive than its comparably sized Random Walk pore scale counterpart, the largest simulation heretofore examined consisted of just under 2,000,000 cells and measured a mere 2.6 mm x 2.6 mm x 1.7 mm. Using a desktop computer with a 2Ghz AMD processor and 4Gb of RAM, diffusion through the model for a simulated 2 hours took 20-40 hours of computation. If it were possible to simulate the entire 10 mm x 10 mm x 20 mm dataset at this resolution, the computation times required to complete the 75 hour lab diffusion experiment would be very large indeed. With the further addition of more components and reactions to the model, it can be deduced that images at this resolution are entirely too large for simulations using this model. Realizing these challenges, three courses of action are being pursued. As has been determined in the study above, binning of the dataset can sometimes be undertaken with very little apparent loss of information or change in the transport characteristics of the model. With this in mind, entire large datasets may be binned down to acceptable sizes before modeling more chemically complex models in the future. The decreased physical complexity will allow for computation times to remain relatively low while the chemical complexities are added and examined. Secondly, smaller volumes of the materials can be examined experimentally. With the above datasets, cropping of the original dataset is achievable because of the lack of 70 laboratory-derived concentration profile data for comparison and verification. Due to the 3D nature of the pore structure and subsequent diffusive transport paths, cropping would not be applicable if such a comparison between simulated and experimental results were to be made. Assumptions of zero mass flux along the lateral boundaries would rightly be called into question. With a smaller experimental dataset, the experimentally set boundary conditions could be used in the simulations and direct comparisons could subsequently be made. Thirdly, the proliferation of multiple-processor computing opens the possibility of creating an enhanced diffusion-specialized parallel computing model. Splitting the simulations into manageable pieces to thus be simulated in parallel amongst multiple processors, computational runtimes can be reduced from their linearly processed counterparts. As should be apparent, these three courses of action are not mutually exclusive. Thus, to quickly examine large, high resolution, chemically complex systems, all three of these actions may be implemented at the same time. Once computational constraints have been addressed, future considerations may consder increasing the complexity of the simulations. In addition to the pure diffusion driven transport, chemical processes may be of importance in contaminated systems. Along with regular surface-contaminant sorption interactions, the implications of other precipitation- dissolution reactions may become prevalent at the small scale. With its dependence on porosity, the diffusive transport paths may be altered significantly as dilution enlarges and precipitation contracts the pore spaces and throats. If widespread precipitation occurs, the transport pathways could even be shut down altogether as the pores all become clogged. As the models become more reliable, future resources may also be put towards examinations of in situ investigation techniques for dataset collection. As in Chitale, 2005, in situ porosity 71 imaging techniques may allow for the models to simulate transport through relatively undisturbed material. This should always be a consideration when dealing with materials from stressed environments, even more so when the stresses are high as they potentially are in deep geologic repositories. With in situ testing, advances in the characterization of potential waste burial and containment sites could be made in the future. 72 Chapter 5: Summary and Conclusions Diffusion, as a means of contaminant transport, is often overlooked or relegated to play a minor role in groundwater systems. As such, it has not been as thoroughly studied as its advective system counterparts. Studies of diffusion in high risk environments, such as geologic repositories considered for the emplacement of nuclear waste, have been relatively recent. In these higher risk environments, contaminant transport must be kept to a minimum as there must be full containment for a perceptively long time. To ensure contaminant containment over these long timelines, studies have begun to examine the small scale processes involved in diffusion-dominated hydraulic settings. With transport being limited to diffusion, the potential for contaminant emission is diminished. To examine the processes of diffusive transport on a local (sub-cm) scale, this study attempts to simulate diffusion at the μm-mm scale using material parameters obtained through high-resolution imaging. With the transport module of the REV scale reactive transport code MIN3P, diffusion was simulated through a small section of the provided material as the two datasets of varying radiometric resolution were subjected to changing voxel resolutions. Varying the radiometric and pixel resolutions, the examination of the potential effects of differing imaging techniques and their effectiveness at certain measurement scales was possible. Unfortunately, the lack of useable target concentration data prevents validation of the simulations to real world results and limits the discussion to comparisons within the simulated results. In the case of 3D diffusion, the use of a binary dataset to imitate pore structure of a pore scale model gives some insight into the potential movement of a contaminant through the material. It can be seen that, due to wandering connections of areas of open porosity, the diffusion of the contaminant through the material is noticeably tortuous and three- 73 dimensional. With decreasing porosity within a material, this reliance on pore connectivity will logically increase, as the amount of open pores gets fewer and farther between. Moving away from the binary dataset, this reliance on connections of pore spaces noticeably decreases. With the increase in bit depth brought about by binning or the simulations involving the 16 bit porosity map, the three-dimensional nature of the diffusion paths becomes less observable as a more unified diffusion front begins to emerge. This smearing of the contaminant into areas of low porosity exists because of the increased bit depth of the 16 bit and binned binary datasets. This increased number of unique porosity values creates a more normal distribution of porosity values and in doing so, decreases the dichotomy between the open pore spaces and the closed matrix. Within the 16 bit dataset, the porosity distribution is already fairly normally distributed and thus, subsequent binning has relatively little effect on the simulations. This lack of effect can also be seen in the constancy of the mass influx into the 16 bit system as there appears to be little change in the connectivity of the pores during the binning process. The largest effect of the binning was seen with the low porosity binary set. This dataset, with its large areas of low porosity matrix material, has the transport occurring through tenuously connected pathways of open pore spaces. When adjacent voxels are averaged arithmetically during the binning process, the increase in bit depth comes at the expense of the low porosity areas. More transport can thus occur through those areas and the tenuous connections between the open pores widens. As the bit depth increases and the porosity values move towards a more normal distribution, mass transfer into the system is also found to increase. In the most extreme case of the low porosity binary Crop 2 being binned by 16 times, the increase in mass transfer is by a factor of 5.26. As the areas of low porosity that originally stop the contaminant from entering the system are binned into areas 74 of medium porosity, the mass can enter the system through the entire input face instead of entering only through areas of open pores. Increased mass transfer into the system can also be seen when the model uses arithmetic averaging instead of harmonic averaging to calculate the mass flux between cells. Though showing similar mass transfer results to the partially binned dataset, the diffusion into the binned dataset remains logical since it is diffusing into an area of non-zero porosity. Contrarily, diffusion into the unbinned dataset using the arithmetically averages mass fluxes transfers mass from the contaminant source and open pores into areas of very low porosity. This increases the amount of mass into the system, but much of the extra mass enters the matrix where it should not belong. Due to this unsuitable transport, harmonic averaging should be used. This is common practice for the solution of groundwater flow problems; however, most solute transport models use arithmetic averaging for calculation of the effective diffusion coefficients. Erroneous mass transfer within the system can also be attributed to errors caused by numerical diffusion. As noted in section 3.3, diffusion numbers in excess of 1 introduce errors caused by the transport of the solute outside of the constraints of the numerical cell in question. When the diffusion number is on the order of 1, the errors are relatively small and occur only at the low concentrations existing at the diffusion front. But as the diffusion number gets sufficiently large (150), errors of 5% occurred even among the bulk of the diffusing mass. Therefore, if computational constraints allow, the diffusion number should be kept below a maximum of 10 to avoid errors that are greater than 1%. Looking at the mass transfer through the systems, there appears a mixture of two forms of transport. Moving from the binary to the binned and 16 bit models, there is a shift from behaving like a pore-scale model with fingers of diffusing material to a model that 75 relies on a more average diffusion front that moves through the system. The former model gives insight into how connectivity of the pores can play a role in the tortuous movement of the solute, but it is highly variable and dependent on sample and instrument resolutions. With the variations seen during the adjustments to the binary dataset, it is difficult to have confidence that the simulations adequately describe the solute transport. The 16 bit porosity model fits into our REV scale definition and produces consistent diffusion profiles due to its more even distribution of porosity. If available, experimentally derived diffusion profiles would be desired to present a form to which the different models could be validated. For future studies, the 16 bit porosity map method of image processing should be further tested and verified against experimentally derived concentrations. As it was found to give the most consistent results through all of its related datasets, it appears to be quite susceptible to data manipulation. As is apparent by the large amounts of data, some sort of data manipulation is needed in order to be able to simulate a sizeable volume of material. If validation attempts find that the 16 bit dataset is inadequately less accurate at predicting diffusion concentration profiles within the material, differing binning methods should be examined. While binning using the arithmetical averaging of adjacent cells appears to give acceptable representations of the original datasets in some respects, it does not appear to work well on others. Notably, the binning of the low porosity dataset allowed for an unreasonable loss of areas of lower porosity. This loss of the low porosity matrix was the cause of increasingly excessive transport through the material as the binning increased. Binning using harmonic averaging would address this problem, but would instead favour the low porosity areas and the pores would instead close up, stopping all transport altogether. Further study into the effects of different binning schema can perhaps find a middle ground. Adhering to the normal averaging formulae, using a geometric mean would perhaps give a 76 desirable result that is the best of both worlds. Looking elsewhere, image resampling techniques aim to manipulate image resolutions while preserving image quality and may provide a means to better binning. Without doubt, in order to move forward with the modeling of diffusion systems at such small scales, computational restraint will be needed to keep the model achievable on today's technology. With advancements in parallel computing, these restrictions may well be loosened. Until then, averaging techniques may allow the laboratory experiments to be conducted such that they can be reliably modeled. 77 References Al-Raoush, R.I., and Willson, C.S. 2005. Extraction of physically realistic pore network properties from three-dimensional synchrotron X-ray microtomography images of unconsolidated porous media systems. J. Contam. Hydrol., 77: 67-89. Altman, S.J., Uchida, M. Tidwell, V.C. Boney, C.M. and Chambers, B.P. 2004. Use of X-ray absorption imaging to examine heterogeneous diffusion in fractured crystalline rocks. J. Contam. Hydrol., 69:1-26. Bear, J. 1979. Hydraulics of Groundwater. McGraw-Hill, New York. Boving, T.B. and Grathwohl, P. 2001. Tracer diffusion coefficients in sedimentary rocks: correlation to porosity and hydraulic conductivity J. Contam. Hydrol., 53: 85–100 Chitale, D.V. 2005. Borehole imaging in reservoir characterization: implementation of a standard interpretation workflow for the clastic- and carbonate reservoirs. SPWLA 46th Annual Logging Symposium, 2005 Finnish Energy Industries, 2006. Nuclear Waste Management in Finland Greiner, A., Schreiber, W. Brix, G. and Kinzelbach, W. 1997. Magnetic resonance imaging of paramagnetic tracers in porous media: Quantification of flow and transport parameters. Water Resour. Res., 33:1461-1473. Gelhar, L.W., Welty, C. and Rehfeldt, K.R. 1992. A Critical Review of Data on Field-Scale Dispersion in Aquifers. Water Resour. Res,. 28:1955-1974. Harvey, C. and Gorelick, S. M. 2000. Rate-Limited mass transfer or macrodispersion: Which dominates plume evolution at the Macrodispersion Experiment (MADE) site? Water Resour. Res,. 36:637-650. Intera, 2006. OPG’s Deep Geologic Repository for Low and Intermediate Level Waste. Geoscientific Site Characterization Plan. INTERA Engineering Ltd. INTERA 05-220-1, OPG 00216-REP-03902-00002-R00. International Energy Agency, 2008. Key World Energy Statistics. http://www.iea.org/publications/free_new_desc.asp?pubs_ID=1199 McDonald, M.G., and Harbaugh, A.W. 1983. A modular three-dimensional finite-difference ground-water flow model. Open-File Report 83-875. U.S. Geological Survey. http://pubs.er.usgs.gov/usgspubs/ofr/ofr83875. Mayer, K. U., Frind, E.O. and Blowes, D.W. 2002. Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water Resour. Res., 38, 78 Mazurek, M. 2004. NWMO Background Paper 6-12, 99 pp. Millington, R.J., and J.M. Quirk. 1960. Transport in porous media. In F.A. Van Beren, et al. (ed.) Trans. 7th Int. Congr. Soil Sci. Vol. 1. 97–106. Nakashima, Y. 2000. The use of X-ray CT to measure diffusion coefficients of heavy ions in water-saturated porous media. Eng. Geol., 56: 11–17. Nakashima, Y. 2003. Diffusivity measurement of heavy ions in Wyoming montmorillonite gels by X-ray computed tomography J. Contam. Hydrol., 61: 147-156. Nakashima, Y., Nakano, T. Nakamura, K. Uesugi, K. Tsuchiyama A. and Ikeda, S. 2004. Three-dimensional diffusion of non-sorbing species in porous sandstone: computer simulation based on X-ray microtomography using synchrotron radiation J. Contam. Hydrol., 74: 253–264. Nakashima, Y., Kamiya, S., and Nakano, T. 2008. Diffusion ellipsoids of anisotropic porous rocks calculated by X-ray computed tomography-based random walk simulations. Water Resour. Res., 44, W12435 Neville, C.J. 2006. Analytical Solutions for Solute Transport with One-Dimensional Flow: Brick Sources in an Infinite Aquifer. Personal communication. Ogata, A and Banks, R.B. A solution of the differential equation of longitudinal dispersion in porous media, U.S. Geol. Surv. Prof. Pap., 411-A, 1961. Oswald, S., Kinzelbach, W. Greiner, A. and Brix, G. 1997. Observation of flow and transport processes in artificial porous media via magnetic resonance imaging in three dimensions. Geoderma, 80:417-429. Patriarche, D., Michelot, J.-L. Ledoux, E. and Savoye, S. 2004. Diffusion as the main process for mass transport in very low water content argillites: 1. Chloride as a natural tracer for mass transport—Diffusion coefficient and concentration measurements in interstitial water. Water Resour. Res., 40:(1) Polak, P., Grader, A.S. Wallach R. and Nativ, R. 2003. Chemical diffusion between a fracture and the surrounding matrix: Measurement by computed tomography and modeling. Water Resour. Res., 39:(4), 1106. Ruiz de Argandona, V.G. Rodriguez Rey, A. Celorio, C. Suarez del Rio, L.M. Calleja L. and Llavona, J. 1999. Characterization by Computed X-Ray Tomography of the Evolution of the Pore Structure of a Dolomite Rock During Freeze-Thaw Cyclic Tests. Phys. Chem. Earth (A), 24:(7), 633-637. Schwartz, F.W. and Zhang, H. 2003. Fundamentals of Ground Water. John Wiley and Sons, Inc. New York, NY, USA 79 Steefel, C.L., MacQuarrie, K.T.B. 1996 Approaches to modeling of reactive transport in porous media. Reactive Transport in Porous Media 34:83-129 Sudicky, E.A., Wadsworth, T.D., Kool, J.B., and Huyakorn, P.S.. 1988. PATCH3D, Three- dimensional analytical solution for transport in a finite thickness aquifer with first type rectangular patch source. Tidwell, V.C., Meigs, L.C. Christian-Frear, T. Boney, C.M. 2000. Effects of spatially heterogeneous porosity on matrix diffusion as investigated by X-ray absorption imaging. J. Contam. Hydrol. 42: 285– 302. Van Geet, M. Lagrou, D. and Swennen, R. 2003. Porosity measurements of sedimentary rocks by means of microfocus X-ray computed tomography (muCT). Geol. Soc., London, Spec. Pub., 215, 51-60. Voss, C.I. and Souza, W.R. 1987. Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone. Water Resour. Res., 23:(10), 1851-1866. Wersin, P. Van Loon, L.R. Soler, J.M. Yllera, A. Eikenberg, J. Gimmi, T. Hernan, P. and Boisson, J.-Y. 2004. Long-term diffusion experiment at Mont Terri: first results from field and laboratory data. Appl. Clay Sc., 26: 123–135. Yllera, A., Hernandez, A. Mingarro, M. Quejido, A. Sedano, L.A. Soler, J.M. Samper, J. Molinero, J. Barcala, J.M. Martın, P.L. Fernandez, M. Wersin, P. Rivas P. and Hernan, P. 2004. DI-B experiment: planning, design and performance of an in situ diffusion experiment in the Opalinus Clay formation. Appl. Clay Sc., 26: 181– 196. Zhang, H., Schwartz, F.W., Wood, W.W., Garabedian, S.P., and LeBlanc, D.R. 1998. Simulation of variable-density flow and transport of reactive and nonreactive solutes during a tracer test at Cape Cod, Massachusetts. Water Resour. Res., 34:(1), 67-82. 80 Appendices A.1 Additional 2D Verification Given the difficulties in finding a robust analytical model code for which to compare diffusion profiles obtained from a 2-dimensional patch source, numerous numerical simulations were run with varying parameters to determine where divergences between certain solutions exist. For the case where the diffusion coefficient equals 1.157d-10 m2/s, subsequent simulations were conducted to determine if boundary effects were present. Table A.1 lists the simulations that were run. Simulation number Column length (m) Spatial increment Δx (m) Column width (m) Spatial increment Δy (m) Source Width Y (m) Diffusion coefficient (m2/s) 48 0.10 0.40 0.03 49 0.20 0.40 0.03 51 0.10 0.80 0.03 55 0.20 0.0025 0.80 0.01 0.03 D0 = 1.157d-10 Table A.1. Simulation parameters used for 2D diffusion while testing for boundary effects. While the spatial increments were kept constant, the length and width of the column were increased. These increases extend the boundaries such that the diffusing Iodide is not influenced by reaching the sides or end of the column. The difference these changing extents have on the diffusion profiles can be seen in Figure A.1. From the data, the concentration reaches the background minimum of 1.0d-10 at a distance along the column of 0.168 m. Thus there should be observable boundary effects for the simulations where the column is only 0.1 m long. 81 -6.E-06 -5.E-06 -4.E-06 -3.E-06 -2.E-06 -1.E-06 0.E+00 0.E+00 1.E-02 2.E-02 3.E-02 4.E-02 5.E-02 6.E-02 7.E-02 8.E-02 9.E-02 1.E-01 Distance along Column Centreline (m) C on ce nt ra tio n D iff er en ce fr om S im ul at io n N 48 (m ol /l) N48 N49 N51 N55 Figure A.1. Concentration difference from original numerical simulation (N48). Diffusion of Iodide along the column centerline. As seen in Figure A.1, there is no observable difference between simulations of similar column lengths, but there is a difference between the shorter and longer column simulations. This difference can be accounted for by the boundary effect that occurs as the Iodide reaches the end of the shorter column. At a distance of 0.02 m along the column centreline, the difference between the simulations of different column lengths is on the order of 1x10-6 mol/l. Comparatively, the difference between the simulated and analytical diffusion profiles at this distance is on the order of 2x10-2 mol/l, indicating that there are fundamental differences between the results that cannot be attributed to boundary effects. 82 This difference can also be seen in the 2D Figure A.2 below. The diffusion pattern is fairly similar between the simulation types, except near the source zone. Figure A.2. Concentration map of Iodide diffusion form a patch source at x=0 for a) MIN3P numerical model, and b) Patch3d analytical model. Near the source zone, the solute appears to be spreading in the y-dimension parallel to the boundary only in the numerical simulation. In the analytical solution, the solute appears to curl around the top and bottom of the source zone and travel towards x=0 at an angle that is not parallel to the y axis. This loss of mass in the –x direction, could account for the lower 83 than simulated analytical concentrations, and is not consistent with the underlying governing equations. 84 A.2 Data Pre-processing The original datasets were supplied from UNB as a set of images stacked along the z- axis, with x and y-coordinates specified for each image. To facilitate input into the MIN3P model, the data needs to be pre-processed to match the code’s data structure. At the same time, the original datasets can be cropped, binned, and transformed to select sub-regions and to adjust the resolution of the data set. Pre-processing was performed using the image manipulation program imageJ (http://rsbweb.nih.gov/ij/) and a set of tailored Fortran 90 programs, developed specifically for this project. For the sake of simplicity, the following example input file is used for each of the created F90 programs. The parameters used below are as follows: nx, ny, nz: number of cells in the x, y, z dimension scale: scaling factor applied to normalize the image to porosity (either 255 or 65536) slice1: the number of the slice from the original dataset that is to be nearest the source change: the distance in metres between slices (Δz) nout: the number of outputs provided by MIN3P nbin: the amount of binning 85 Example Input nx ny nz scale slice1 change 144 144 96 255 300 0.180E-04 nout 1 nbin 2 A.2.1 Cropping Cropping of the stack of images in the x-y plane is completed using the built-in cropping tool within the image manipulation program imageJ (http://rsbweb.nih.gov/ij/). ImageJ allows for cropping of multiple images at a time. The required starting pixels and size of crop is all that is required. As pictured in Figure IV.1, Crop 1 is located 20,20 pixels from the top left and is 144x144 pixels. Crop 2 is located 215,215 from the top left corner of the original dataset and is also 144x144 pixels. Cropping of the stack in the z direction is completed by choosing the desired number of slices to enter into the transformation program (Appendix A.2.2.) A.2.2 Transformation into MIN3P Format In order to correctly read the input data, the datasets are transformed from their original x-y image slices into a single file of the format x,y,z,n. This program reads in the stack dimensions (nx,ny,nz), the amount of scaling to be applied to each voxel (scale), the location of the first slice to be read (slice), and the distance in metres between the slices (change). Also inputted are each x-y slice named prefix####.txt where the prefix is specified on running the program and the number increases sequentially through the slices. The program outputs the transformed file as output.out in the format x,y,z,n without the required three line header. The header needs to be added manually after the transformation. The 86 program also calculates and outputs the average value of n per slice for each slice in average.avg. PROGRAM generalized INTEGER :: nx INTEGER :: ny INTEGER :: nz INTEGER :: nn INTEGER :: nc INTEGER :: nt INTEGER :: xlow INTEGER :: xhigh REAL :: xtemp REAL :: scale = 1 REAL :: xcng REAL :: ycng REAL :: zcng = 0 REAL :: chng INTEGER :: ivol = 1 INTEGER :: jvol = 1 INTEGER :: kvol = 1 INTEGER :: lvol = 1 INTEGER :: tvol = 1 INTEGER :: iz = 1 INTEGER :: itxt = 10 INTEGER :: iout = 15 INTEGER :: iavg = 12 INTEGER :: iinp = 11 INTEGER :: fstslc = 0 character*72 :: prefix character*4 :: suffix character*1 :: cdummy INTEGER :: l_prfx INTEGER :: l_sufx REAL, PARAMETER :: pi = 3.141592 REAL, allocatable :: por(:) REAL, allocatable :: por_slice(:) REAL, allocatable :: rowt(:,:) REAL, allocatable :: row(:,:) REAL, allocatable :: av(:) REAL, allocatable :: xnum(:) REAL, allocatable :: ynum(:) REAL, allocatable :: znum(:) write(*,*)'enter the prefix of the files to be concatenated...' c getting prefix name call getstrqq(prefix) l_prfx=index(prefix,' ')-1 if(l_prfx.eq.-1.or.l_prfx.gt.72)then l_prfx=72 endif 87 c read dimensions from input file open(iinp,file=prefix(:l_prfx) & //'.inp',status='unknown',form='formatted') read(iinp,*) cdummy read(iinp,*) nx,ny,nz,scale,fstslc,chng c read(iinp,*) nx,ny,nz,scale,chng c allocating parameters nn = nx * ny nt = nn * nz allocate (por(nt)) allocate (por_slice(nn)) allocate (rowt(1,nx)) allocate (row(nx,1)) allocate (av(nz)) allocate (xnum(nt)) allocate (ynum(nt)) allocate (znum(nt)) c start z loop do iz = 1,nz lvol = fstslc + iz-1 c lvol = 300 + iz-1 ycng = 0 100 format(I4.4) write (suffix,100) lvol l_sufx = 4 open(itxt,file=prefix(:l_prfx)//suffix(:l_sufx) & //'.txt',status='unknown',form='formatted') c transposing file kvol = 1 do jvol = 1,ny read(itxt,*) (row(ivol,1),ivol=1,nx) rowt=transpose(row) xcng = 0 do ivol=1,nx if ((rowt(1,ivol))==0) then por_slice(kvol)=(0.1000000E-19) else por_slice(kvol)=((rowt(1,ivol)))/scale end if por(tvol)=por_slice(kvol) xnum(tvol)=xcng ynum(tvol)=ycng 88 znum(tvol)=zcng xcng = xcng + chng kvol = kvol + 1 tvol = tvol + 1 end do ycng = ycng + chng end do zcng = zcng + chng c averaging of slices av(iz) = sum(por_slice) / size(por_slice) close (itxt) end do c end z do loop open(iavg,file='average.avg',status='unknown',form='formatted') write(iavg,'(/a)') 'average:' write(iavg,'(72a/)')('-',i=1,72) write(iavg,'(/a)') 'the average for the specified slice & in percent(%) is:' do iz = 1,nz write(iavg,*) 100 * av(iz) end do close (iavg) open(iout,file='output.out',status='unknown',form='formatted') do ivol = 1,nt write(iout,200) xnum(ivol),ynum(ivol),znum(ivol),por(ivol) end do 200 format(100E15.7) close (iout) END PROGRAM generalized A.2.3 Binning Like cropping, binning of the stack of images in the x-y plane is also completed using the built-in imageJ plugin filter. Binning in these planes was completed as consecutive 2x binnings. After binning in the x-y plane, the dataset is transformed as explained above. 89 Binning of the stack in the z direction is completed after the transformation using the following program. This program reads in the stack dimensions (nx,ny,nz) as well as the desired binning degree (nbin) from the input file. Also inputted is the transformed dataset file prefix.por with the three line header from the transformation step. This program outputs a dataset file named prefix.binned in the same x,y,z,n format as was inputted. PROGRAM binXtimes INTEGER :: nx INTEGER :: ny INTEGER :: nz INTEGER :: zbin INTEGER :: nnew INTEGER :: nn INTEGER :: nt INTEGER :: nout INTEGER :: need INTEGER :: nbin = 2 INTEGER :: ivol = 1 INTEGER :: tvol = 1 INTEGER :: iout = 1 INTEGER :: ix = 1 INTEGER :: iz = 1 INTEGER :: ibinned = 12 INTEGER :: iinp = 11 INTEGER :: ipor = 10 character*72 :: prefix character*1 :: cdummy INTEGER :: l_prfx INTEGER :: l_sufx REAL, allocatable :: binz(:) REAL, allocatable :: por(:) REAL, allocatable :: por_sum(:) REAL, allocatable :: xnum(:) REAL, allocatable :: ynum(:) REAL, allocatable :: znum(:) write(*,*)'enter the prefix of the files to be binned...' c getting prefix name call getstrqq(prefix) l_prfx=index(prefix,' ')-1 if(l_prfx.eq.-1.or.l_prfx.gt.72)then l_prfx=72 endif c read dimensions from input file 90 open(iinp,file=prefix(:l_prfx) & //'.inp',status='unknown',form='formatted') read(iinp,*) cdummy read(iinp,*) nx,ny,nz read(iinp,*) cdummy read(iinp,*) cdummy read(iinp,*) cdummy read(iinp,*) nbin c allocate sizes to arrays nx = (nx / nbin) ny = (ny / nbin) zbin = (nz / nbin) nn = nx * ny nt = nn * nz nnew = nn * zbin allocate (binz(nnew)) allocate (por(nt)) allocate (por_sum(nn)) allocate (xnum(nt)) allocate (ynum(nt)) allocate (znum(nt)) c read in porosity array open(ipor,file=prefix(:l_prfx) & //'.por',status='unknown',form='formatted') c read(ipor,*) cdummy c read(ipor,*) cdummy c read(ipor,*) cdummy do ivol = 1,nt read(ipor,*) xnum(ivol),ynum(ivol),znum(ivol),por(ivol) end do close (ipor) c loop to get average per slice do z = 0,(zbin-1) c initialize por_sum array do ivol = 1,nn por_sum(ivol) = 0.0 end do do b = 0,(nbin-1) ivol = 1 do y = 0,(ny-1) do x = 1,nx need = (nn * ((z * nbin) + b)) + (ny * y) + x 91 por_sum(ivol) = por_sum(ivol) + por(need) if (b .eq. (nbin-1)) then binz(tvol) = por_sum(ivol) / nbin tvol = tvol + 1 end if ivol = ivol + 1 end do end do end do end do c write average output file for each average output open(ibinned,file=prefix(:l_prfx) & //'.binned',status='unknown',form='formatted') write(ibinned,'(/a)') 'normalized average:' write(ibinned,'(72a/)')('-',i=1,72) write(ibinned,'(/a)') 'the normalized concentration is:' do iz = 1,(nnew) write(ibinned,200) xnum(iz),ynum(iz),znum(iz),binz(iz) end do 200 format(100E15.7) close (ibinned) END PROGRAM binXtimes 92 A.3 Modifications to MIN3P To accommodate the needs of the datasets provided for the study, numerous modifications had to be made to MIN3P. Those modifications were completed to allow for variable archie`s law exponents (Equation 1.5, Section 1.2), to replace the arithmetic averaging scheme with a harmonic scheme for diffusion coefficients (Equations 2.5, 2.5, Section 2.3.4), and to allow for the input of a porosity field. A.3.1 Archie`s Law Exponent To account for a variable Archie`s law exponent, the parameter had to be read from the input file and subsequently used for the determination of the diffusion coefficient. First, the initcprt subroutine was modified to initialize and then search for and read the exponent from the input file as follows: initcprt.f parameter (r86400 = 8.64d4, r3 = 3.0d0, r4 = 4.0d0) archie = r4/r3 c tortuosity corrections - new definition subsection = 'tortuosity correction' call findstrg(subsection,itmp,found_subsection) if (found_subsection) then read(itmp,*,err=999,end=999) subsection if (subsection.eq.'millington') then tortuosity_corr = .true. elseif (subsection.eq.'archie') then tortuosity_corr = .true. if (found_subsection) then read(itmp,*,err=999,end=999) archie end if elseif (subsection.eq.'no correction') then tortuosity_corr = .false. end if end if 93 Once within the program, the exponent was passed to the diffcoff subroutine as the variable arch to calculate the diffusion coefficient. diffcoff.f real*8 function diffcoff(diff_p,sat_p,por,tortuosity_corr,arch) implicit real*8 (a-h,o-z) parameter (r3 = 3.0d0, r4 = 4.0d0, r10 = 10.0d0) logical tortuosity_corr if (tortuosity_corr) then diffcoff = diff_p * sat_p**(r10/r3) * por**(arch) else diffcoff = diff_p * sat_p * por end if return end The diffcoff subroutine is called in a number of subsequent subroutines and they also had to be modified to call the new arch variable. These included infcrt_a, infcrt_g, jacbrt, and mbalrt. A.3.2 Harmonic Averaging of Diffusion Coefficients The harmonic averaging of the diffusion coefficients is computed within the subroutine infcrt_a_d. This was repurposed from the original subroutine infcvs which dealt with the harmonic averaging of permeabilities. It was modified to read and work with the diffusion coefficients. 94 infcrt_a_d.f c ---------------------------------------------------------------------- c subroutine infcvs tweaked for diffusion c ----------------- c c compute influence coefficients (variably saturated flow) c c c /| ij, K_ij c / | / c / |/ c | / c K_i | /| K_j c *--------|-*-------------* c i d_i,ij | \| d_j,ij j c | \ c | / \ c | / A_ij c c c written by: Uli Mayer - June 12, 96 c c last modified: Scot Ellis - January 8, 07 c c ---------------------------------------------------------------------- subroutine infcrt_a_d use parm use gen use phys implicit real*8 (a-h,o-z) parameter (rhalf = 0.5d0, r1 = 1.0d0) sa = r1 !fully saturated c initialize pointer to current control volume ivol = 0 c loop over control volumes do ivz = 1,nvz ! number of control volumes in z do ivy = 1,nvy ! number of control volumes in y do ivx = 1,nvx ! number of control volumes in x c pointer to current control volume ivol = ivol+1 jtemp = iavs(ivol) c assign conductivities for current control volume 95 izn = mpropvs(ivol) !retrieve material properties por_i = pornew(ivol) !assign porosity diff_i = diffcoff(diff_a,sa,por_i,tortuosity_corr, & archie) !assign diff coef c permeability update due to porosity changes c if (update_permeability) then c cxx_i = perm_fac(ivol)*cxx_i c cyy_i = perm_fac(ivol)*cyy_i c czz_i = perm_fac(ivol)*czz_i c end if jtemp = jtemp+1 ! skip diagonal c pointers to previous colums in x,y and z ivxp = ivx-1 ivxn = ivx+1 ivyp = ivy-1 ivyn = ivy+1 ivzp = ivz-1 ivzn = ivz+1 c assign interfacial distances of current control volume if (half_cells) then !half_cells on boundary if (nvx.gt.1) then !in x-direction if (ivx.eq.1.or.ivx.eq.nvx) then !boundary delx_i = delx(ivx) elseif (ivx.gt.0.and.ivx.lt.nvx) then !interior delx_i = rhalf*delx(ivx) end if end if if (nvy.gt.1) then !in y-direction if (ivy.eq.1.or.ivy.eq.nvy) then !boundary dely_i = dely(ivy) elseif (ivy.gt.0.and.ivy.lt.nvy) then !interior dely_i = rhalf*dely(ivy) end if end if if (nvz.gt.1) then !in z-direction if (ivz.eq.1.or.ivz.eq.nvz) then !boundary delz_i = delz(ivz) elseif (ivz.gt.0.and.ivz.lt.nvz) then !interior delz_i = rhalf*delz(ivz) end if end if else !full cells on boundary if (nvx.gt.1) then !in x-direction 96 delx_i = rhalf*delx(ivx) end if if (nvy.gt.1) then !in y-direction dely_i = rhalf*dely(ivy) end if if (nvz.gt.1) then !in z-direction delz_i = rhalf*delz(ivz) end if end if !(half_cells) c calculate influence coefficients in x-direction if (nvx.gt.1) then !connections in x-direction if (ivxp.gt.0) then !left connection (2) areaf = dely(ivy)*delz(ivz) if (half_cells) then !half cells on boundary if (ivxp.eq.1) then delx_j = delx(ivxp) elseif (ivxp.gt.1) then delx_j = rhalf*delx(ivxp) end if else !full cells on boundary delx_j = rhalf*delx(ivxp) end if !(half_cells) jvol = javs(jtemp) jzn = mpropvs(jvol) por_j = pornew(jvol) diff_j = diffcoff(diff_a,sa,por_j,tortuosity_corr, & archie) c permeability update due to porosity changes c if (update_permeability) then c cxx_j = perm_fac(jvol)*cxx_j c end if diff_ij = diff_i*diff_j cinfrt_da(jtemp) = diff_ij*areaf/(diff_i*delx_j+ & diff_j*delx_i) jtemp = jtemp+1 end if if (ivxn.le.nvx) then !right connection (3) areaf = dely(ivy)*delz(ivz) if (half_cells) then !half cells on boundary if (ivxn.eq.nvx) then delx_j = delx(ivxn) elseif (ivxn.lt.nvx) then delx_j = rhalf*delx(ivxn) end if else !full cells on boundary delx_j = rhalf*delx(ivxn) 97 end if !(half_cells) jvol = javs(jtemp) jzn = mpropvs(jvol) por_j = pornew(jvol) diff_j = diffcoff(diff_a,sa,por_j,tortuosity_corr, & archie) c permeability update due to porosity changes c if (update_permeability) then c cxx_j = perm_fac(jvol)*cxx_j c end if diff_ij = diff_i*diff_j cinfrt_da(jtemp) = diff_ij*areaf/(diff_i*delx_j+ & diff_j*delx_i) jtemp = jtemp+1 end if end if !connections in x-direction c calculate influence coefficients in y-direction if (nvy.gt.1) then !connections in y-direction if (ivyp.gt.0) then !front connection (4) areaf = delx(ivx)*delz(ivz) if (half_cells) then !half cells on boundary if (ivyp.eq.1) then dely_j = dely(ivyp) elseif (ivyp.gt.1) then dely_j = rhalf*dely(ivyp) end if else !full cells on boundary dely_j = rhalf*dely(ivyp) end if !(half_cells) jvol = javs(jtemp) jzn = mpropvs(jvol) por_j = pornew(jvol) diff_j = diffcoff(diff_a,sa,por_j,tortuosity_corr, & archie) c permeability update due to porosity changes c if (update_permeability) then c cyy_j = perm_fac(jvol)*cyy_j c end if diff_ij = diff_i*diff_j cinfrt_da(jtemp) = diff_ij*areaf/(diff_i*dely_j+ & diff_j*dely_i) jtemp = jtemp+1 end if if (ivyn.le.nvy) then !back connection (5) 98 areaf = delx(ivx)*delz(ivz) if (half_cells) then !half cells on boundary if (ivyn.eq.nvy) then dely_j = dely(ivyn) elseif (ivyn.lt.nvy) then dely_j = rhalf*dely(ivyn) end if else !full cells on boundary dely_j = rhalf*dely(ivyn) end if !(half_cells) jvol = javs(jtemp) jzn = mpropvs(jvol) por_j = pornew(jvol) diff_j = diffcoff(diff_a,sa,por_j,tortuosity_corr, & archie) c permeability update due to porosity changes c if (update_permeability) then c cyy_j = perm_fac(jvol)*cyy_j c end if diff_ij = diff_i*diff_j cinfrt_da(jtemp) = diff_ij*areaf/(diff_i*dely_j+ & diff_j*dely_i) jtemp = jtemp+1 end if endif !connections in y-direction c calculate influence coefficients in z-direction if (nvz.gt.1) then !connections in z-direction if (ivzp.gt.0) then !bottom connection (6) areaf = delx(ivx)*dely(ivy) if (half_cells) then !half cells on boundary if (ivzp.eq.1) then delz_j = delz(ivzp) elseif (ivzp.gt.1) then delz_j = rhalf*delz(ivzp) end if else !full cells on boundary delz_j = rhalf*delz(ivzp) end if !(half_cells) jvol = javs(jtemp) jzn = mpropvs(jvol) por_j = pornew(jvol) diff_j = diffcoff(diff_a,sa,por_j,tortuosity_corr, & archie) c permeability update due to porosity changes c if (update_permeability) then c czz_j = perm_fac(jvol)*czz_j 99 c end if diff_ij = diff_i*diff_j cinfrt_da(jtemp) = diff_ij*areaf/(diff_i*delz_j+ & diff_j*delz_i) jtemp = jtemp+1 end if if (ivzn.le.nvz) then !top connection (7) areaf = delx(ivx)*dely(ivy) if (half_cells) then !half cells on boundary if (ivzn.eq.nvz) then delz_j = delz(ivzn) elseif (ivzn.lt.nvz) then delz_j = rhalf*delz(ivzn) end if else !full cells on boundary delz_j = rhalf*delz(ivzn) end if !(half_cells) jvol = javs(jtemp) jzn = mpropvs(jvol) por_j = pornew(jvol) diff_j = diffcoff(diff_a,sa,por_j,tortuosity_corr, & archie) c permeability update due to porosity changes c if (update_permeability) then c czz_j = perm_fac(jvol)*czz_j c end if diff_ij = diff_i*diff_j cinfrt_da(jtemp) = diff_ij*areaf/(diff_i*delz_j+ & diff_j*delz_i) jtemp = jtemp+1 end if end if !connections in z-direction end do ! number of increments in x end do ! number of increments in y end do ! number of increments in z cdbg c do irow=1,nn c istart = iavs(irow) c iend = iavs(irow+1)-1 c write(igen,'(6(a,i5,a,f10.3,1x))') c & ('cinfvs(',i1,') = ',cinfvs(i1),i1=istart,iend) c end do c stop cdbg return end 100 A.3.3 Allowing Input of Porosity Field To have the model simulate diffusion through the provided porosity fields, these fields had to be inputted and processed by MIN3P. For this, an array (pornew) was created in the place of the original porosity parameter. The porosity file was then opened using opngfls, read into the program and assigned using initpppm. opngfls ipor = 26 c porosity distribution open(ipor,file=prefix(:l_prfx)//'.por',status='unknown', & form='formatted') initpppm character*1 cdummy parameter (r1 = 1.0d0, tiny = 1.0d-10) c set default porosity_field = .false. if (varsat_flow.or.reactive_transport) then c read porosity field from file subsection = 'read porosity field from file' call findstrg(subsection,itmp,found_subsection) if (found_subsection) then porosity_field = .true. c read porosities read(ipor,*,err=999,end=999) cdummy 101 read(ipor,*,err=999,end=999) cdummy read(ipor,*,err=999,end=999) cdummy do ivol = 1,nn read(ipor,*,err=999,end=999) rdummy, rdummy, rdummy, & por_init(ivol) end do close (ipor) end if c read porosity from input file if (.not.porosity_field) then read(icnv,*,err=999,end=999) por end if if (porosity_field) then pornew(ivol) = por_init(ivol) else pornew(ivol) = por por_init(ivol) = por end if nntemp = nntemp+1 mpropvs(ivol) = izn !allocate properties Since the original porosity parameter was now an array, some loops had to be constructed within subroutines batreac, initbcrt, initicrt, and transbcrt. Batreac is shown as a sample. batreac do ivol = 1,nn !loop over control volumes if (porosity_field) then swc = 1.0d0 sac = 0.0d0 else call rtrvpprm(swc,sac,pornew(ivol),section_header) endif 102 c compute initial condition call gcreact(ccnew,ccold,cxc,gamma_l(1),gamma_l(nc+1), & cgc,swc,sac,pornew(ivol),igen,ilog,idbg, & tec_header,prefix,l_prfx,zone_name,l_zone_name) end do !loop over control volumes c compute initial condition do ivol = 1,nn !loop over control volumes call gcreact(ccnew,ccold,cxc,gamma_l(1),gamma_l(nc+1), & cgc,swc,sac,pornew(ivol),igen,ilog,idbg, & tec_header,prefix,l_prfx,zone_name,l_zone_name) end do !loop over control volumes A.4 Dataset Post-processing Once the simulation is complete, the following program reads in the concentration data from the simulation output from multiple output times and calculates and outputs the average concentration per slice for each slice as prefix_##.avg where ## corresponds to the output time number. From the general input file above, the program reads in the stack dimensions (nx,ny,nz) as well as the number of output time files to average (nout). PROGRAM profile_norm INTEGER :: nx INTEGER :: ny INTEGER :: nz INTEGER :: nn INTEGER :: nc INTEGER :: nt INTEGER :: nout INTEGER :: need INTEGER :: ivol = 1 INTEGER :: tvol = 1 INTEGER :: inorm = 1 INTEGER :: iout = 1 INTEGER :: ix = 1 INTEGER :: iz = 1 INTEGER :: iavg = 12 INTEGER :: iinp = 11 INTEGER :: igst = 13 103 INTEGER :: ipor = 10 INTEGER :: inor = 14 character*72 :: prefix character*4 :: suffix character*1 :: cdummy c character*1 :: borp = 'b' INTEGER :: l_prfx INTEGER :: l_sufx REAL, allocatable :: av(:) REAL, allocatable :: avconc(:) REAL, allocatable :: conc(:) REAL, allocatable :: conc_need(:) REAL, allocatable :: avpor(:) REAL, allocatable :: por(:) REAL, allocatable :: por_need(:) REAL, allocatable :: por_conc(:) REAL, allocatable :: norm(:) REAL, allocatable :: xnum(:) REAL, allocatable :: ynum(:) REAL, allocatable :: znum(:) write(*,*)'enter the prefix of the files to be averaged...' c getting prefix name call getstrqq(prefix) l_prfx=index(prefix,' ')-1 if(l_prfx.eq.-1.or.l_prfx.gt.72)then l_prfx=72 endif c read dimensions from input file open(iinp,file=prefix(:l_prfx) & //'.inp',status='unknown',form='formatted') read(iinp,*) cdummy read(iinp,*) nx,ny,nz read(iinp,*) cdummy c read(iinp,*) nout,borp read(iinp,*) nout c allocate sizes to arrays nn = nx * ny nt = nn * nz allocate (av(nz)) allocate (avconc(nz)) allocate (conc(nt)) allocate (conc_need(nn)) allocate (avpor(nz)) allocate (por(nt)) allocate (por_need(nn)) allocate (por_conc(nn)) allocate (norm(nt)) allocate (xnum(nt)) allocate (ynum(nt)) allocate (znum(nt)) 104 c read in porosity array open(ipor,file=prefix(:l_prfx) & //'.por',status='unknown',form='formatted') read(ipor,*) cdummy read(ipor,*) cdummy read(ipor,*) cdummy do ivol = 1,nt read(ipor,*) rdummy, rdummy, rdummy, por(ivol) end do close (ipor) c read in outputted concentration array 100 format (I1) c loop for each output file do iout = 1,nout write (suffix,100) iout l_sufx = 1 open(igst,file=prefix(:l_prfx)//'_'//suffix(:l_sufx) & //'.gst',status='unknown',form='formatted') read(igst,*) cdummy read(igst,*) cdummy read(igst,*) cdummy do ivol = 1,nt read(igst,*) xnum(ivol),ynum(ivol),znum(ivol),conc(ivol) end do close (igst) iz = 1 c loop to get average per slice do z = 0,(nz-1) ivol = 1 do y = 0,(ny-1) do x = 1,nx need = (nn * z) + (ny * y) + x por_need(ivol) = por(need) conc_need(ivol) = conc(need) por_conc(ivol) = (conc_need(ivol) * por_need(ivol)) c if (borp .EQ. 'p') then c norm(inorm) = (conc(inorm) / por(inorm)) c inorm = inorm + 1 c end if ivol = ivol + 1 end do end do 105 106 avpor(iz) = sum(por_need) / size(por_need) avconc(iz) = sum(por_conc) / size(por_conc) c avconc(iz) = sum(conc_need) / size(conc_need) av(iz) = avconc(iz) / avpor(iz) iz = iz + 1 end do !! write (*,*) 100*av(iz) c write average output file for each average output open(iavg,file=prefix(:l_prfx)//'_'//suffix(:l_sufx) & //'.avg',status='unknown',form='formatted') write(iavg,'(/a)') 'Maximum value is:' write(iavg,*) maxval(conc) write(iavg,'(72a/)')('-',i=1,72) write(iavg,'(/a)') 'normalized average:' write(iavg,'(72a/)')('-',i=1,72) write(iavg,'(/a)') 'the normalized average for the specified & slice is:' do iz = 1,(nz) write(iavg,*) av(iz) end do close (iavg) c write full output of normalized concentration data if porosity input c if (borp .EQ. 'p') then c open(inor,file=prefix(:l_prfx)//'_'//suffix(:l_sufx) c & //'.nor',status='unknown',form='formatted') c write(inor,'(/a)') 'normalized average:' c write(inor,'(72a/)')('-',i=1,72) c write(inor,'(/a)') 'the normalized concentration is:' c c do iz = 1,(nz) c c write(inor,200) xnum(iz),ynum(iz),znum(iz),norm(iz) c c end do c200 format(100E15.7) c close (inor) c end if end do END PROGRAM profile_norm
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical modeling techniques for assessing three-dimensional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical modeling techniques for assessing three-dimensional diffusion processes in heterogeneous rock… Ellis, Scotland Russell 2010
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Numerical modeling techniques for assessing three-dimensional diffusion processes in heterogeneous rock samples on the sub-mm scale |
Creator |
Ellis, Scotland Russell |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | The burial of high level nuclear wastes in geologic repositories requires careful consideration of the long-term hydrogeological and geochemical stability of the receiving formations. Care must be taken due to the high environmental sensitivity of the waste material and its potential long-term effect on groundwater. Studies into the host rocks' natural ability to minimize contaminant migration are a matter of high priority in the planning for long-term storage of high level radioactive waste. Focusing on porous sedimentary rock, this study aims to examine numerical techniques used to analyze and interpret experimental data that characterize the distribution of porosity in geologic samples at the µm to mm scale. Because repositories are hosted in natural fine-grained rock formations, transport in the vicinity of these repositories is diffusion-controlled and believed to be affected substantially by heterogeneities at these scales. Data available for the analyses consist of non-destructive, high-resolution measurements of porosity obtained using new developments in X-ray imaging. Advances in computing technology make it possible to numerically analyze the expected patterns of sub-mm-scale diffusive transport for these large experimental data sets. The modeling analyses examine 3D diffusive transport in heterogeneous rock samples and evaluate the effect of data resolution and image processing techniques on the connectivity of the transport pathways. The simulation results provide insight into small-scale diffusive transport of solutes, and guide the needs for dataset resolution and handling for these types of investigations. With increased availability of experimental results, further modeling studies could be conducted. These studies would aim at developing a link between simulation results and observed data to further develop the transport theory for contaminant migration on this scale. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-11-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0052917 |
URI | http://hdl.handle.net/2429/30233 |
Degree |
Master of Applied Science - MASc |
Program |
Geological Engineering |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_spring_ellis_scotland.pdf [ 1.25MB ]
- Metadata
- JSON: 24-1.0052917.json
- JSON-LD: 24-1.0052917-ld.json
- RDF/XML (Pretty): 24-1.0052917-rdf.xml
- RDF/JSON: 24-1.0052917-rdf.json
- Turtle: 24-1.0052917-turtle.txt
- N-Triples: 24-1.0052917-rdf-ntriples.txt
- Original Record: 24-1.0052917-source.json
- Full Text
- 24-1.0052917-fulltext.txt
- Citation
- 24-1.0052917.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}]}"
data-media="{[{embed.selectedMedia}]}"
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-0052917/manifest