TRANSPORT OF SORBING SOLUTES IN FRACTURED MEDIA A NUMERICAL AND EXPERIMENTAL ANALYSIS OF DISPERSION AND RETARDATION by CHRISTOPH WELS B.Sc. Trent University, Peterborough, 1985 M.Sc. Trent University, Peterborough, 1989 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES Department of Geological Sciences We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA APRIL 1995 © Christoph Wels in presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of G^OLOGICDL The University of British Columbia Vancouver, Canada Date DE-6 (2/88) X^fQ^I°lC SC(£/UC££ Abstract The primary pathways for contaminant transport in lowpermeability fractured rock are likely to be through a network of hydraulically-connected fractures in the rock formation. Sorption of contaminants to the fracture walls may significantly retard their transport. The influence of sorption on solute transport is analyzed using both numerical and laboratory migration experiments. A random walk model is developed to simulate solute transport in a parallel plate fracture, assuming that fast, reversible, and linear sorption occurs in a small zone adjacent to the fracture wall. It is shown that a sorbing solute experiences greater longitudinal spreading than does a conservative solute. The magnitude of this enhanced dispersion reaches a maximum in the range of fluid velocities characteristic of Taylor dispersion. At the network scale, transport is simulated by tracking particles in discrete fracture networks, assuming that within each fracture, retardation varies in proportion to the product of a surface distribution coefficient and the specific surface area of the fracture. The results suggest that retardation at the plume scale is a non-uniform and anisotropic process. Different segments of the plume, or equivalently different breakthrough fractions at a downstream boundary, are retarded to a different degree. The degree to which various ii breakthrough fractions are retarded varies as a function of the orientation of the mean hydraulic gradient relative to the orientation of the fracture sets. This variation can be described in the form of a retardation ellipse. To test the surface retardation model used in the numerical analyses, the sorbing radionuclides strontium and uranium were injected in smooth-walled fractures in granite and steel, respectively. The retardation of uranium was inversely proportional to fracture aperture, providing qualitative support for the definition of the surface retardation factor. The influence of fracture aperture on the retardation of strontium was much greater than that predicted by the surface retardation factor. Strontium retardation was approximately an order of magnitude greater in a narrow fracture (b=450jim, Ra=45) compared to that in a wide fracture (b=780jim, Ra=3.5) . It is hypothesized that hysteresis in sorption, in conjunction with limited transverse mixing across the fracture, caused the apparent increase in sorption strength with a decrease in fracture aperture. iii Table of Contents Page Abstract ii Table of Contents iv List of Tables viii List of Figures ix Acknowledgements xii Preface xiii Chapter One: INTRODUCTION 1 Chapter Two: THE INFLUENCE OF SURFACE SORPTION ON SOLUTE DISPERSION IN PARALLEL-PLATE FRACTURES 6 2.1 Introduction 6 2.2 Random Walk Model 12 2.3 Dispersion of a Non-sorbing Solute 20 2.4 Dispersion of a Sorbing Solute 24 2.5 Implications for Nonuniform Retardation 33 2.6 Implications for Solute Transport in Rough-Walled iv Fractures 2.7 36 Implications for Modeling Solute Transport at the Network Scale 38 2.8 Conclusions 40 2.9 References 43 Chapter Three: THE INFLUENCE OF SPECIFIC SURFACE AREA ON SOLUTE TRANSPORT IN FRACTURES - AN EXPERIMENTAL ANALYSIS 52 3.1 Introduction 52 3.2 Materials and Methods 60 3.2.1 Migration Experiments 60 3.2.2 Static Sorption Experiments 64 3.2.3 Interpretative Models 66 3.3 Results and Discussion 3.3.1 71 Influence of Fracture Aperture on Retardation of Uranium 3.3.2 71 Static Sorption Experiments of Strontium on Granite 74 3.3.3 Analysis of Tritium Breakthrough Curves 3.3.4 Influence of Non-Ideal Sorption on Dispersion of Strontium 3.3.5 82 Influence of Fracture Aperture on Retardation of Strontium 3.3.6 . 80 89 Implications for Reactive Transport in Fractured Media 98 v 3.4 Conclusions 100 3.5 Notation 103 3.6 Acknowledgements 104 3.7 References 105 Chapter Four: RETARDATION OF SORBING SOLUTES IN FRACTURED MEDIA 124 4.1 Abstract 124 4.2 Introduction 125 4.3 Retardation Models 132 4.4 Numerical Migration Experiments: Design and Choice 4.5 of Parameters 142 Results and Discussion 151 4.5.1 Two Orthogonal Fracture Sets with a Common Aperture Distribution 4.5.2 A Continuum Model for Nonuniform Retardation? 4.5.3 155 Influence of Fracture Density on Retardation 4.5.4 161 Two Orthogonal Fracture Sets with Differing Mean Aperture 4.5.5 165 Influence of Variability in Fracture Orientation on Retardation 4.5.6 151 Limitations of this Study 174 177 4.6 Conclusions 181 4.7 Acknowledgements 183 vi 4.8 References Chapter Five: 184 CONCLUSIONS AND RECOMMENDATIONS 201 5.1 Concluding Remarks 201 5.2 Summary of Conclusions 20 6 vii List of Tables Table Page 3.1. Experimental conditions in migration experiments . . 108 3.2. Chemical composition of tracer solutions 109 3.3. Summary of surface distribution coefficients . . . . 110 3.4. Summary of hydraulic parameters determined for migration experiments Ill 4.1. Summary of fracture network input parameters . . . . viii 187 List of Figures Figure Page 2.1. Conceptualization of reactive solute transport in a parallel-plate fracture 45 2.2. The fracture Peclet number as a function of mean fluid velocity 46 2.3. Simulated breakthrough curves of a conservative solute 47 2.4. Simulated breakthrough curves of sorbing solutes . . 48 2.5. The ratio Deff/DL as a function of the mean fluid velocity 49 2.6. The dispersion factor cf determined numerically for a range of retardation factors 50 2.7. The ratio Deff/DL and retardation factors Rref as a function of transport distance 51 3.1a. Experimental set-up for continuous injection . . . 3.1b. Experimental set-up for pulse injection 112 113 3.2. Autoradiographs of the uranium plume on the fracture surface 114 3.3. The concentration profile of removable uranium on the fracture surfaces 115 3.4. Experimental results of batch sorption of strontium on crushed granite 116 3.5. Sorption of strontium to three granite coupons as a ix function of time 117 3.6. Desorption of strontium from three granite coupons as a function of time 118 3.7. Experimental breakthrough of tritium 119 3.8. Experimental breakthrough of strontium and model simulations for experiment 3 120 3.9. Comparison of experimental breakthrough of strontium in experiment 4 and 5 121 3.10. Experimental breakthrough of strontium and model type curves for experiment 4 122 3.11. Experimental breakthrough of strontium and model type curves for experiment 5 123 4.1. Flow domain, boundary conditions and the source zone used in the simulations 188 4.2. Single realization of a fracture network representing the reference geometry 189 4.3. Transport behavior in the dense network of scenario 1 190 4.4a. Spatial particle distributions for a single realization of scenario 1 with o equal to 0.5 for the conservative solute plume 191 4.4b/c. Spatial particle distributions for a single realization of scenario 1 with a equal to 0.5 for the sorbing solute plume 192 4.5. Histograms of the particle velocities x 193 4.6. Transport behavior in the sparse network of scenario 2 194 4.7. Flow distribution in single realizations of the dense network of scenario 1 and the sparse network of scenario 2 for G equal to 0.5 195 4.8. Transport behavior in scenario 3 and 4 196 4.9. Retardation factors plotted in polar coordinates for a equal to 0.0 and o equal to 0.5 . 197 4.10. Example realization of a fracture network in scenario 5 and scenario 6 198 4.11. Retardation ellipses for scenario 5 and 6 for o equal to 0.0 199 4.12. Retardation ellipses for scenario 5 and 6 for o equal to 0.5 200 xi Acknowledgements Above all, I wish to thank my thesis advisor Leslie Smith for his invaluable advice and support throughout this study. I am very grateful for his patient guidance in every aspect of this thesis, from the early stages of discussing research ideas to the final editing of each thesis chapter. I also thank Roger Beckie and Tom Brown for their help as members of my thesis committee. Special thanks go to Roger Beckie for his advice on modeling and the world of numbers. This thesis has benefitted greatly from the interaction with my colleagues and friends in the hydrogeology group. Among many others I like to thank Tom Clemo, Bob Parney, Dan Walker, and Brian Berkowitz for the numerous discussions which helped focussing my ideas. I am indebted to Tom Clemo for making his computer code DISCRETE available to me. Special thanks also go to Chuck Vandergraaf for his advice and logistical support during my stay at the Whiteshell Nuclear Laboratories. Finally, I would like to thank my wife Rocio E. Lopez who provided wonderful support and encouragement throughout this study. This research was financially supported by the National Science and Engineering Research Council of Canada and the University of British Columbia's Graduate Fellowship Program. xii Preface In accordance with University guidelines, and with the agreement of my Supervisory Committee, this thesis consists of a series of three research papers, each with a distinct focus, and each forming a chapter of the thesis. The first two of these research papers (Chapters Two and Three) represent manuscripts to be submitted to a scientific journal in the near future. The third research paper (Chapter Four) has been published recently (Wels and Smith, 1994). The published research paper is co-authored by Dr. Leslie Smith, thesis supervisor, who has provided advice, support, and editorial guidance throughout the research. I am grateful to my supervisor and co-author for his continuing advice and assistance. All of the research, analysis, and interpretation presented in the three papers was performed by Christoph Wels, thesis author, in accordance with the guidelines of the University of British Columbia. Publications: Wels, C. and L. Smith, Retardation of sorbing solutes in fractured media, Water Resourc. Res., 30(9), 2547-2563, 1994. Acknowledged: xiii Chapter One: INTRODUCTION Low-permeability geologic formations such as crystalline rock or clayey tills have been proposed as host formations for the isolation of radioactive waste and other hazardous substances. It is reasoned that the low hydraulic conductivity of the host material minimizes the movement of groundwater which may carry dissolved contaminants away from the waste site. These "tight" formations may, however, contain a network of interconnected fractures which could compromise the effectiveness of this natural barrier. In a low-permeability rock, open fractures provide much less resistance to flow than does the rock matrix. As a result groundwater flows primarily in these fractures, carrying with it the dissolved contaminants. The rate of advective transport in the network of fractures may be reduced significantly by chemical interaction of the contaminants with the rock matrix and/or minerals lining the fracture. Clearly, a comprehensive understanding of reactive solute transport in fractured media is fundamental to assessing the risks of hazardous waste disposal in low-permeability geologic formations. In this thesis we study the influence of sorption to the fracture walls on solute transport in fractured media. We 1 focus on the chemical process of sorption for two reasons. First, adsorption/desorption between the solid phase (rock) and the aqueous phase (solution) is known to have a strong influence on the transport of many hazardous substances, in particular trace metals and radionuclides. Second, the influence of sorption on solute transport has been studied in great detail for the case of porous media. It has been shown that many sorbing solutes exhibit the same movement through the porous medium as a non-reactive solute, but with a travel velocity reduced by a (constant) retardation factor. It is unclear, however, whether and how the concept of retardation can be applied to solute transport in fractured media. The concept of retardation is based on experimental data which show that partitioning of many solutes between the solid and aqueous phase is a fast, reversible process which can be described by a linear isotherm. With the assumption of linear equilibrium sorption the transport equation for a sorbing solute in a homogeneous porous medium is given by V-(JPVC) _ V-(CV) R R dC 8t .-. ' l where C is the solute concentration in solution, V is the average linear velocity of the groundwater and D is the dispersion tensor. R is known as the retardation factor and is defined as 2 R = 1 + ^KD (2) where 0 and pb are the porosity and the bulk density of the porous medium, respectively. KD is the bulk distribution coefficient which expresses the extent of sorption relative to the mass of the porous medium available for sorption. It is seen from equation (1) that the effect of sorption is to decrease both the advective flux and the dispersive flux by a constant retardation factor, R. The retardation concept can be extended to transport of sorbing solutes in fractured media. In a fractured medium, sorption may occur during diffusive transport in the unfractured rock matrix (bulk sorption) and/or during advective-dispersive transport in the fractures (surface sorption). For the case of bulk sorption within the rock matrix the retardation model is equivalent to that of a porous medium (equation (2)). For the case of surface sorption to the fracture walls, however, the retardation model is better expressed in terms of a surface retardation factor Ra Ra = 1 + ^Ka (3) where SA/V is the surface area-to-volume ratio of the fracture (or network of fractures) . Ka is defined as the surface distribution coefficient which expresses the extent of sorption relative to the surface area of the rock available for sorption. The objective of this thesis is to evaluate the application of the retardation concept to solute transport in fractured media. In Chapter Two we evaluate the conditions under which the advection-dispersion equation (1) is a valid description of transport of a sorbing solute in a parallel plate fracture. The process of surface sorption is modeled explicitly using a random walk model which assumes that sorption occurs only in the immediate vicinity of the fracture walls. Numerical migration experiments demonstrate that a one-dimensional form of the advection-dispersion equation (1), in conjunction with the surface retardation model (3), accurately describes solute transport under the assumption that there is sufficient transverse mixing of the solute across the fracture aperture. It will also be shown that the process of surface sorption results in enhanced longitudinal spreading of the solute plume and nonuniform retardation for "short" fracture lengths. In Chapter Three we test the surface retardation model (3) for transport of the sorbing solutes strontium and uranium in constant-aperture fractures machined from granite and steel, respectively. The radionuclides were injected in smooth-walled fractures that differed only in fracture aperture. Results of these laboratory migration experiments demonstrate the strong influence of the specific surface area (=SA/V) of a fracture 4 on transport of a sorbing solute. It will be shown that the surface retardation model (3) may significantly underestimate solute retardation if the sorption process is not fully reversible. In Chapter Four we evaluate the application of the retardation concept to model solute transport at the network scale. The movement of a sorbing solute is simulated in discrete fracture networks assuming the retardation model (3) holds at the scale of a single fracture. Numerical migration experiments show that at the plume scale, retardation is both nonuniform and anisotropic. It will be shown that the degree of nonuniformity and anisotropy in the retardation factor is largely controlled by the fracture geometry and the orientation of the fracture network in the flow field. Our results suggest that a field-scale continuum model based on the conventional advection-dispersion equation (1), using a uniform retardation factor (3), may be conceptually incorrect, even for a densely fractured medium. In the concluding chapter recommendations are given for future research on transport of sorbing solutes in fractured media. Chapter Five closes with a summary of the main conclusions which were drawn in the present study. 5 Chapter Two: THE INFLUENCE OF SURFACE SORPTION ON SOLUTE DISPERSION IN PARALLEL-PLATE FRACTURES 2.1 Introduction Fracture network models commonly represent each fracture in the network as parallel plates with constant aperture (e.g., Endo et al., 1984/ Long and Witherspoon, 1985; Clemo, 1994). Solute transport is modeled explicitly in each fracture using either an advection-dispersion equation (e.g., Sudicky and McLaren, 1992) or assuming transport by advection alone (e.g., Schwartz et al., 1983/ Andersson and Dverstorp, 1987). The former approach assumes that solute dispersion in each fracture is "Fickian". The latter approach assumes that: (i) dispersion is negligible relative to advection in each fracture, and/or (ii) dispersion at the single fracture scale is negligible relative to the dispersion at the network scale (e.g., Hull et al., 1987). Network models can be extended to account for sorption to the fracture walls through the introduction of a surface retardation factor (e.g., Sudicky and McLaren, 1992/ Dverstorp et al., 1992/ Wels and Smith, 1994). The use of a constant retardation factor for a given 6 fracture implicitly assumes that there is no interaction between surface sorption and solute dispersion in a single fracture. Dispersion of a conservative solute in a parallel-sided fracture is governed to a large extent by the fracture geometry. In an open fracture of constant aperture b, solute disperses under the combined influence of the parabolic velocity profile between the fracture walls and molecular diffusion, a process known as Taylor dispersion (Taylor, 1953; 1954) . Beyond an entrance length Lcrit required to eliminate concentration variations perpendicular to flow, dispersion is "Fickian" and transport can be approximated by the advectiondispersion equation with a longitudinal dispersion coefficient DL given as (Aris, 1956) D u2b2 =D + L m 210Dm 11) ' K where u is the mean fluid velocity in the fracture and Dm is the molecular diffusion coefficient. Note that hydrodynamic dispersion in a parallel-plate fracture is caused by diffusion in the longitudinal direction (second term in (1)) as well as Taylor dispersion (third term in (1) ) . The entrance length Lcrlt for this asymptotic dispersion coefficient to apply is given as (Kessler and Hunt, 1994) 7 _ 6 ub2 _ Lerit (2) ~~^^7 Hull et al. (1987) used the same relationship with 0.25 as proportionality constant (instead of 6/n2) to determine the applicability of the Fickian dispersion model in constantaperture fractures. Assuming a fast, linear, and reversible sorption onto the fracture walls solute transport in a parallel-sided fracture can be described by the one-dimensional form of the advectiondispersion equation (e.g., Tang et al. , 1981) *.f->i-»i where C is the concentration of the solute in solution, DL is the longitudinal dispersion coefficient (equation (1)), and Ra is the surface retardation factor defined as Ra = l+ ^ (4) Ka is known as the surface distribution coefficient which can be estimated from static sorption experiments (e.g., Vandergraaf et al, 1988) . The use of a constant retardation factor in the transport equation (3) assumes that retardation 8 of the solute plume is uniform, that is, each segment of a plume is retarded to the same degree. Equation (3) further assumes that there is no interaction between surface sorption and solute dispersion. In other words, the magnitude of DL is assumed to be independent of the sorption strength Ka of the sorbing solute. Little is known about the interaction between surface sorption and solute dispersion in parallel-sided fractures. Most laboratory migration experiments which have used machined, constant-aperture fractures in rock have indicated that sorbing solutes exhibit significantly more dispersion than conservative solutes (e.g., Vandergraaf et al., 198 8; Holtta et al., 1991; Wels, 1995). Possible reasons for this enhanced dispersion include (i) spatial variability in sorption sites/potential, (ii) non-equilibrium sorption, and/or (iii) other non-ideal sorption effects such as nonlinear or partially irreversible sorption (e.g., Wels, 1995). However, enhanced dispersion is not limited to the case of non-ideally sorbing solutes. Kessler and Hunt (1994) used analytical expressions to show that (ideal) sorption in a porous surface coating of a parallel-sided fracture may also enhance longitudinal dispersion of the solute plume. Note that enhanced dispersion due to sorption results in nonuniform retardation at the plume scale, regardless of which mechanism is causing it. In this chapter we focus on the interaction between sorption 9 to the fracture walls and solute dispersion in a parallelplate fracture. Although a parallel-plate fracture does not reproduce all features of natural fractures (e.g., channeling or surface roughness), the parallel-plate conceptualization has proven very useful in the study of non-reactive solute transport in fractures, and it is thus a legitimate starting point for the reactive case. For simplicity, we assume an impermeable host rock so that diffusion into and sorption within the matrix can be neglected. Furthermore, surface sorption is assumed to be a fast, linear and reversible process. Our analysis focuses on the localized nature of the sorption process (at the fracture walls) and its influence on solute dispersion. This study is based on the working hypothesis that solute has to move into the vicinity of the fracture wall in order to participate in the sorption process. The above hypothesis would limit the application of the transport equation (3) in describing transport of a sorbing solute (Ra>l) in a parallel-sided fracture. In the case of a conservative solute (Ra=l) the one-dimensional transport equation may be applied once the solute has traveled far enough into the fracture to establish transverse homogenization of the solute by molecular diffusion. In the case of a solute, however, that sorbs to the fracture walls a second condition must be invoked. It is necessary to assume that the reactive solute mixes completely and instantaneously across the fracture aperture. Note that this assumption of 10 "instantaneous transverse mixing" for a sorbing solute is much more limiting than that of "transverse homogenization" for a conservative solute. Transverse homogenization gradually occurs as the solute moves along the fracture. A transport distance greater than the entrance length Lcrlt defined in equation (2) is sufficient to homogenize the solute front of a conservative solute. In contrast, the applicability of instantaneous transverse mixing is controlled by the rate of transverse diffusion relative to advective transport along the fracture and is independent of the transport distance. The objective of this study is to evaluate the conditions under which the advection-dispersion equation (3) is a valid description of reactive transport in a parallel-plate fracture, and to evaluate limitations in its application. We present numerical migration experiments which suggest that, in comparison to a non-sorbing solute, a sorbing solute requires a greater entrance length into the fracture before transverse homogenization is established. In addition, surface sorption is shown to result in significant enhanced dispersion under transport conditions which favour advective transport along the fracture over transverse diffusion across the fracture aperture. An effective dispersion coefficient is derived to account for this enhanced dispersion due to surface sorption. 11 2.2 Random Walk Model In a parallel-plate fracture (Fig. 2.1a) solutes are transported by two mechanisms: (i) advection by the fluid which travels along the fracture at different velocities according to the Poiseuille velocity profile and (ii) molecular diffusion. For a fracture aspect ratio greater than approximately 12, the transport problem is essentially twodimensional which can be written as (e.g., Hull et al., 1987) dt m{ ax2 dy2' 2 { K } a 2 ' dx where a is the half-aperture (=b/2) of the fracture (Fig. 2.1b). The process of fast, linear, and reversible sorption to the fracture walls is described by the following boundary condition Cs = CKa \y=a (6) where C is the dissolved solute concentration (mass per unit volume) in proximity of the surface and Cs is the surface solute concentration (mass sorbed per unit area) on the fracture wall. We employ a random walk method to simulate the movement of conservative and sorbing solutes in a parallel-plate fracture 12 represented by the above transport equations (5) and (6). This method is well-suited to study dispersion phenomena in a fracture because the processes of molecular diffusion (probabilistic), advection (deterministic), and sorption (probabilistic) are modeled explicitly. The transport of a solute by advection and diffusion is modeled by moving a large number of particles according to a random walk with steps of the form (e.g., Kinzelbach, 1987) xp(t+At) = xp(t) + Z1^2D^At + v(yp)At (7a) and yp(t+At) =yp(t) + Z2J2D^M (7b) where xp and yp are the longitudinal and transverse position of a particle, respectively, At is the time step and Zlf Z2 are random normal variates with mean zero and variance one. The local particle velocity v(yp) is determined from the parabolic velocity profile of the fluid which is given as v(yp) = | U (i-4) <8> In the limit of large numbers of particles and small time steps, the frequency distribution obtained from the random 13 walk (7a,b) satisfies the conservative transport equation (5). This methodology has been extended to account for surface sorption onto the fracture wall. A "sorption zone" is defined within a small distance sw of the fracture wall (Fig. 2.1c) in which surface sorption is assumed to take place. The use of a (relatively small) sorption zone near the fracture wall is consistent with the working hypothesis in this chapter that surface sorption influences only those solute particles in close proximity of the sorbing fracture wall. In our conceptualization sorption occurs (strictly speaking) only at the fracture surface (see boundary condition (6)). From a theoretical point-of-view a solute would have to either hit the surface directly for chemi-sorption to occur or at least move within the thickness of the electric double layer for surface complexation or ion exchange to occur. In either case, the distance from the surface within which a solute particle may sorb, is probably not greater than several tens of nanometers. Such a discretization of the sorption problem was not feasible in the particle transport code. A larger sorption width on the order of the distance a particle may travel by diffusion in a few time steps was considered acceptable. The use of a larger sorption width tends to underestimate the degree of enhanced dispersion due to surface sorption. The partitioning of solute between solid phase and solution is expressed by the linear isotherm 14 [g] Ka = -&A [C] O) where [q] and [C] are the molar concentrations of the solute sorbed onto the surface and in solution, respectively. For use in a particle-based random walk, we need to express equation (9) in terms of the number of particles (representing mass of solute) assigned to each phase. The equivalent partitioning of particles between the surface and solution of the sorption zone is given by the partitioning coefficient Kp ^n = — P S„ = — (10) N In words, the partitioning of particles is proportional to Ka and the surface area-to-volume ratio (SA/V) of the sorption zone (i.e. l/sw in a two-dimensional description). For a given sorption zone, a greater sorption strength (large Ka) results in a greater number of particles Ns being sorbed relative to particles N present in the sorption zone sw. This partitioning of particles between surface and solution is modeled probabilistically in the random walk method. Every particle is assigned an additional state variable sp(t) indicating whether the particle is in the dissolved state (sp(t)=0) or in the sorbed state (sp(t)=l). Every time step, particles residing within the sorption zone (i.e. yp is within 15 a distance sw from the fracture wall; see Fig. 2.1c) are evaluated for a transition between the dissolved and the sorbed state according to the following transition probabilities Po,i = •£%- ; Pi,0 = j ~ (11) p P These transition probabilities describe equilibrium sorption and are consistent with the more general transition probabilities for a first-order rate controlled sorption reaction given by Valocchi and Quinodoz (1989). The use of the equilibrium transition probabilities (11) assumes that the time scale of sorption is much smaller than the time scale of advection (local equilibrium). For a particle in the dissolved state at time t, the particle is moved according to the transport step (7a,b) before determining the new state variable according to sp(t+At) = 0 for X > p01 1 for X < p 0fl (12) where X is a uniformly distributed random variable in the interval [0,1]. In contrast, for a particle in the sorbed state at time t there is no movement, i.e. the transport step 16 is replaced by a sorption step xp(t+At) =xp(t) } yp(t+At) = yp(t) (13) and the new state variable is then determined according to s (t+At) = 0 for X sp 1 for X > p lf0 (14) In the limit of a large number of particles (Np), a small time step (At) and a small sorption width (sw) , the frequency distribution obtained from the transport step (7a,b) combined with the sorption step (11)-(14) simulates the two-dimensional transport with surface sorption onto the fracture wall expressed by equations (5)-(6). For each time step the particle is first moved by diffusion to its new coordinates. Based on the old and new y-coordinate, an average particle velocity in the x-direction is calculated and the particle is moved according to this average fluid velocity. Particles may intersect the fracture wall and center-line only in the diffusive step. Both boundaries represent no-flow boundaries which reflect particles back into the model domain. For particles that reflect off the fracture wall or the centerline, the average particle velocity is given 17 as the average velocity for the path from the old position to the boundary and the path from the boundary to the new position. Note that the time step At is chosen small enough so that a particle can not leave the sorption zone when reflected off the fracture wall. This way all particles which "hit" the fracture wall are automatically evaluated for sorption using equation (11). Particles are released into the fracture across the upgradient boundary of the fracture as a Dirac pulse. A flowweighted (flux) injection is obtained by releasing the particles in proportion to the flow rate at any given transverse position. The transport parameters supplied to the model are four physical parameters (u,Dm,b,L) and one chemical parameter (Ka) . In addition the algorithm requires the use of three model parameters (Np,At,sw) . In general, the larger Np and the smaller the values of At and sw the smoother the breakthrough curve (BTC) and the more accurate the solution. Sensitivity analyses indicated that 2000 particles and a time step small enough to obtain a ratio of sorption width to diffusion distance (sw/(2DmAt)1/2) greater than four is sufficient to simulate the transport problems discussed here. It was also found that the reactive transport solution is fairly insensitive to the relative width (a/sw) of the sorption zone. As a result of the shape of the parabolic velocity profile (high fluid velocities concentrated close to the center of the fracture) the value of sw may be as large as 18 15 percent of the half aperture without significantly influencing the reactive transport solution. The non-dimensionalized form of the advection-dispersion equation (3) is given as _ dC R 1 d2C dC ( 1 5 ) Pe dX2 dX dT where the dimensionless variables are defined as X= X L Pe = uL m R - (16) ut LRf In our application, L is the transport distance of interest, DL is the longitudinal dispersion coefficient defined in equation (1), and Pe is the Peclet number. For conservative transport (Ra=l) the dimensionless time TR simplifies to T=ut/L and all other non-dimensional variables remain unchanged. For a continuous injection the analytical solution to the non-dimensional advection-dispersion equation (15) is (Ogata and Banks, 1961) X T 1+ T Pe = | eerfc(— rfc( * ))++±e—Pee erfc( erf c (— R* 2 R 2 o 2^T /Pe 2^T /Pe — C 19 ) (17) Because (17) is linear in C the solution for a continuous injection is equivalent to the solution for a pulse injection integrated over time (Kreft and Zuber, 1978). Hence, equation (17) can be fit to the non-dimensionalized cumulative recovery curves obtained with the random walk model. The cumulative recovery curves eliminate much of the stochastic variations present in the instantaneous BTCs, thus simplifying a least squares fitting procedure. Note that only one fitting parameter (Pe) is needed due to the non-dimensionalization of the problem. 2.3 Dispersion of a Non-sorbing Solute In this section we will briefly review the transport of a non-sorbing solute in a parallel-plate fracture which will serve as reference for the discussion of the reactive transport case. Assuming the solute has traveled farther than the entrance length Lcrit required to establish transverse homogenization, transport of a conservative solute ina parallel-plate fracture is described by the advectiondispersion equation (15) with the following fracture Peclet number Pe UL u2b2 £> + 210Dm (18) 20 The numerator describes the magnitude of advection and the denominator describes the magnitude of hydrodynamic dispersion in a parallel-plate fracture. Note that longitudinal dispersion in a parallel-plate fracture consists of longitudinal dispersion due to molecular diffusion and Taylor dispersion. Three transport regimes may be distinguished on the basis of the relative magnitude of the two mechanisms of dispersion accounted for in the denominator of equation (18). In order to delineate these transport regimes we define the transverse Peclet number as Pe t = -5- (19) m When the transverse Peclet number is less than approximately two the transport regime exhibits "diffusion-controlled dispersion". Under these conditions Taylor dispersion accounts for less than two percent of the longitudinal spreading. In contrast, when the transverse Peclet number is greater than 100, longitudinal diffusion accounts for less than two percent of all longitudinal spreading. The latter transport regime is dominated by Taylor dispersion and the effects of longitudinal diffusion may be ignored. For intermediate values of Pet, both mechanisms of dispersion are of similar magnitude and each process needs to be considered explicitly. This transport regime is said to exhibit "mixed dispersion". 21 The fracture Peclet number (18) indicates a complex dependence of conservative solute transport on the fluid velocity u across the three transport regimes. Figure 2.2 shows Pef values as a function of u for b=400jim, Dm=2xl0"9 m2s_1 and several values of transport distance (L=0.02, 0.2 and 2.0m). The vertical dashed lines delineate the three transport regimes. In the transport regime of diffusion-controlled dispersion (u < lxlO-5 ms-1) the fracture Peclet number is proportional to fluid velocity (Fig. 2.2). In the transport regime of Taylor dispersion (u > 5xl0-4 ms-1) the Peclet number is inversely proportional to fluid velocity (Fig. 2.2). The region of mixed dispersion is a transition zone where the fracture Peclet number is comparatively insensitive to changes in u. Note that the highest fracture Peclet numbers are realized in the regime of mixed dispersion (Fig. 2.2). In other words, the non-dimensional breakthrough curves tend to be the steepest for intermediate fluid velocities characteristic of the regime of mixed dispersion. These general observations hold for any transport distance L of interest (Fig. 2.2) provided L is greater than the entrance length Lcrit. These analytical results were confirmed using simulations with our random walk model. Figure 2.3 shows simulated breakthrough curves of a conservative solute for transport conditions representing diffusion-controlled dispersion (circles) , mixed dispersion (triangles) , and Taylor dispersion 22 (squares). The three transport simulations differ only with respect to fluid velocity (covering four orders of magnitude in u ) . All other physical parameters are constant (L=0.2m, b=400|am/ Dm=2xl(T9 m2s-1) . The lines in Figure 2.3 represent the analytical solution (17) using the respective Peclet number calculated from equation (18). As anticipated, Figure 2.3 shows that the transport simulations agree well with the analytical solutions based on transverse homogenization, for the cases of diffusion-controlled dispersion and mixed dispersion. In the Taylor regime, with u=5xl0~3 ms _1 , the analytic approximation of Kessler and Hunt (1994) gives an entrance length, Lcrit, equal to 0.24m. The match in the Taylor regime is also good (for L=0.2 m ) . Thus, our numerical results are consistent with analytic predictions across all three dispersive regimes. Furthermore, they suggest that the entrance length proposed by Kessler and Hunt provides a conservative approximation for transverse homogenization to occur. Figure 2.3 also demonstrates the relative importance of dispersion compared to advection for the three transport regimes. The transport simulation representing mixed dispersion exhibits a steeper breakthrough curve (larger Pef) than the transport simulations representing either Taylor dispersion or "diffusion-controlled dispersion" (Fig. 2.3). In other words, the ratio of dispersive transport to advective transport is the lowest in the regime of mixed dispersion. 23 This simulated transport behaviour is accurately predicted by the fracture Peclet number (18) provided the transport distance allows for transverse homogenization (equation 2 ) . 2.4 Dispersion of a Sorbing Solute In this section we evaluate under which conditions the reactive transport model (3)-(4) is valid in describing transport of a solute that sorbs to the fracture walls in a parallel-plate fracture, and limitations in its application. Two questions deserve special attention. First, under what condition is transverse homogenization of a sorbing solute established? Second, what is the effect of surface sorption on dispersion? Numerical simulations have been carried out to address these questions. Figure 2.4 shows simulated breakthrough curves of two solutes with differing sorption strengths (Ka=2.0xl0~4 and 1.0xl0~2m), representing "weak" retardation (Ra=2) and "strong" retardation (Ra=51), for otherwise identical transport conditions (u=5xl0~3ms-1, L=0.2m, b=400jim, Dm=2xl0"9m2s""1) . The BTC for a conservative solute is shown for comparison (Ra=l). Note that the breakthrough curves are plotted relative to a dimensionless time for the reactive solute TR, that is, the dimensionless time T=ut/L is normalized to the surface retardation factor (TR=T/Ra) . 24 Figure 2.4 illustrates three key effects of surface sorption on solute transport in a planar fracture. First, the greater retardation due to a higher sorption strength is described accurately by the surface retardation factor (4). In other words, all three BTCs shown in Figure 2.4 have the same mean arrival time equal to one in dimensionless time TR. Second, the sorption strength of a solute influences the entrance length required to establish transverse homogenization. In our example simulations (Fig. 2.4) a transport distance of 0.2m is sufficient to establish transverse homogenization of the conservative solute, as indicated by the good fit of the advection-dispersion equation (solid line). In the case of the weakly-sorbing solute, the transport model gives a satisfactory fit, suggesting that transverse homogenization is nearly established (dashed line, Fig. 2.4). In contrast, the same fracture length is not sufficient to establish transverse homogenization of the strongly-sorbing solute, as indicated by the poorer fit of the retardation model (dotted line, Fig. 2.4). The transport distance had to be increased by more than an order of magnitude (L>10m) to obtain a good fit with the advectiondispersion model. These results and additional sensitivity analyses suggest that the increase in the entrance length Lcrlt is of the same magnitude as the surface retardation factor for a given solute. Hence we propose to extend the definition of the entrance length derived by Kessler and Hunt (1994) for the 25 conservative case (equation (2)) to Lent = K « - ^ H £ - (20) Recall that the entrance length indicates the transport distance necessary to eliminate concentration variations across the fracture aperture. The process of surface sorption slows down the opportunity for transverse homogenization because a fraction of the solute mass is sorbed to the surface and can not participate in advective-diffusive transport. With an increase in sorption strength, more mass is sorbed to the surface at any point in time and the rate of transverse mixing by diffusion is further reduced. As a result, a strongly sorbing solute requires a greater transport distance to establish transverse homogenization. Third, Figure 2.4 demonstrates that surface sorption leads to enhanced dispersion, that is, there is a greater spread in the reactive BTC relative to the conservative BTC (Fig. 2.4). This enhanced dispersion is due to the localized nature of surface sorption. In our random walk model solute particles can only participate in the sorption process when they are in the proximity of the fracture wall, that is, when they are present in the sorption zone sw. Those solute particles outside the sorption zone cannot sorb. The localized nature of the sorption reaction thus increases the range of residence 26 times for solute particles travelling in the fracture, and hence, enhances dispersion. Furthermore, the parabolic velocity profile introduces a negative cross-correlation between local fluid velocity and the potential for sorption. There is a higher probability for a slow-moving solute to sorb than for a fast-moving solute because it is closer to the fracture wall, or more likely to be within the sorption zone sw. This bias in sorption further increases the spreading of the solute plume. The extent to which dispersion is enhanced by surface sorption can be quantified by determining an effective dispersion coefficient that is based on the Peclet number obtained from the fit of the advection-dispersion equation to the simulated BTC. Note that this effective dispersion coefficient Deff is always greater than the longitudinal dispersion coefficient DL which describes hydrodynamic dispersion of a conservative solute. In the above example, the effective dispersion coefficient calculated for the weaklysorbing solute is Deff=7.14*1CT5 m2s_1 (Pe=14.0) compared to DL=9.5*10~4 m2s_1 (Pe=105) for a conservative solute. In other words, the effective dispersion coefficient for the weaklysorbing solute is approximately 7.5 times greater than the longitudinal dispersion coefficient (i.e. Deff/DL=7.5) . The effective dispersion coefficient for the strongly-sorbing solute should be even greater. However, for the stronglysorbing solute, a value of Deff could not be estimated in this 27 example because the advection-dispersion model does not provide a good fit to the simulated BTC. The estimation and subsequent application of Deff is only meaningful if the transport distance is sufficiently long to ensure homogenization of the sorbing solute in question (see equation (20)) . The ratio of Deff to DL is a surrogate for the increase in dispersion due to surface sorption beyond the (known) hydrodynamic dispersion in a parallel-plate fracture. Figure 2.5 shows the ratio of Deff to DL for the case of a weaklysorbing solute (Ra=2) over a wide range of fluid velocities. First, consider the solid dots which represent values of Deff/DL determined numerically using the random walk model. The numerical simulations suggest that enhanced dispersion due to surface sorption is insignificant (Deff/DL=1.0) in the transport regime of diffusion-controlled dispersion. In contrast, in a transport regime dominated by Taylor dispersion, longitudinal spreading due to the surface sorption process is much greater than that due to hydrodynamic dispersion alone (Deff/DL»1.0) . In this regime the ratio Deff/DL reaches a plateau value which is nearly independent of fluid velocity (Fig. 2.5). We will show below that this value depends on the magnitude of the retardation factor. For a retardation factor of 2, Deff/DL is approximately equal to 7.8 in the Taylor regime which is close to the value (Deff/DL=7.5) obtained from the simulation shown in Figure 2.4 (for R a =2). At intermediate fluid velocities the 28 ratio Deff/DL gradually increases from unity to the plateau value, indicating a strong dependence of the magnitude of enhanced dispersion on fluid velocity in the transport regime of mixed dispersion (Fig. 2.5) . These numerical results are best understood by re-examining the origin of the enhanced dispersion. Longitudinal spreading of the solute plume increases in the presence of surface sorption because only those particles in the vicinity of the fracture wall may sorb. The negative cross correlation between local fluid velocity and the probability for sorption further enhances solute dispersion. This additional dispersion is absent only if sorption does not depend on the transverse position of the solute. In theory, this is accomplished when the solute particles mix instantaneously across the fracture. In practice, it is sufficient to require that transverse diffusion is much faster than advection. This condition is met in the transport regime of diffusion-controlled dispersion (Pet<2). Here, the fluid velocities are so small that diffusion dominates over advection and instantaneous transverse mixing is approximated. It is important to distinguish the condition of "instantaneous transverse mixing" of a solute from that of "transverse homogenization". Transverse homogenization is gradually achieved as the solute front moves along the fracture. In contrast, the applicability of instantaneous transverse mixing is controlled by the transport regime 29 (Pet<2) and it is independent of the transport distance. For transport conditions that do not approximate instantaneous transverse mixing (Pet>2), enhanced dispersion occurs beyond transport distances sufficient to produce transverse homogenization. Sensitivity analyses suggested that the magnitude of enhanced dispersion is only influenced by the magnitude of Taylor dispersion. This observation is expressed in the following relationship ».* " D. • ct ^ (21) where the proportionality factor cf is defined here as a "dispersion factor". The above formulation only applies for transport distances sufficiently large to allow for transverse homogenization of the sorbing solute. The solid line in Figure 2.5 shows the behaviour of Deff/DL as a function of the mean fluid velocity, using the above definition of the effective dispersion coefficient (21) with cf=7.8. This line fits the simulated data well in all three different transport regimes (Fig. 2.5). Hence the above formulation of an effective dispersion coefficient appears general and applies to a wide range of transport conditions. Note that cf is best determined in the Taylor regime. In this 30 region longitudinal dispersion due to molecular diffusion is negligible, and as a result, the ratio Deff/DL approximates the dispersion factor (Fig. 2.5). The same value for the dispersion factor (cf=7.8) was obtained in simulations with other combinations of b, v and Dm suggesting that the magnitude of the dispersion factor depends on the magnitude of Taylor dispersion (=u2b2/210Dm) and not on the individual magnitude of the physical parameters b, v and Dm. Figure 2.6 shows the dispersion factors cf as a function of surface retardation factors. The data points were obtained from transport simulations using the random walk model. The results shown in Figure 2.6 apply in all transport regimes, provided the transport distance L is greater than the entrance length required to achieve transverse homogenization of the sorbing solute (equation (20)). For conservative transport (Ra=l) the dispersion factor is by definition equal to one. Figure 2.6 shows that a higher retardation factor results in a larger dispersion factor indicating more enhanced dispersion due to surface sorption. For high surface retardation factors, however, the dispersion factor seems to approach an asymptotic value (cf~20) . The increase in the magnitude of enhanced dispersion with increasing sorption strength is consistent with the proposed causes of the enhanced longitudinal spreading. It has been suggested that only a limited fraction of the solute particles (those in immediate vicinity of the fracture wall) participate 31 in surface sorption. With an increase in the sorption strength a greater proportion of the solute mass that is within sw participates in this limited and biased sorption process and the magnitude of enhanced dispersion increases. The dispersion factor approaches a maximimum for strongly-sorbing solutes (high Ra) where the great majority of the solute particles present in the sorption zone sw are sorbed to the surface. Here an incremental change in sorption strength has little effect on the partitioning of the solute mass and the resulting enhanced dispersion. In summary, transport of a sorbing solute in a parallelplate fracture may be modeled using the advection-dispersion equation (3) with the surface retardation factor (4) under conditions that allow for transverse homogenization of the sorbing solute. The dispersion coefficient describing hydrodynamic dispersion needs to be modified, however, to account for additional dispersion due to surface sorption. The definition of the effective dispersion coefficient (equation (21)) is applicable for a wide range of transport conditions since the dispersion factor is only a function of the sorption strength of the solute (Fig. 2.6). The magnitude of enhanced dispersion due to surface sorption increases with sorption strength. The effective dispersion coefficient for a stronglysorbing solute may be up to 20 times greater than the longitudinal dispersion coefficient for a conservative solute. 32 2.5 Implications for Nonuniform Retardation It has been shown that surface sorption results in enhanced longitudinal spreading of the solute plume in a parallel-plate fracture. At the plume scale, this enhanced dispersion due to surface sorption leads to nonuniform retardation, that is, the leading segment of the plume experiences less retardation than the center of the plume and the tail of the plume experiences more retardation. In the following we will demonstrate that the degree of nonuniform retardation depends on the distance that a solute plume travels in a fracture. The influence of transport distance on nonuniform retardation was studied by comparing the breakthrough curves of a conservative solute and a sorbing solute (Ra=5.0), simulated using the random walk model for a range of transport distances (L=0.01-lm). For any given transport distance we computed three retardation factors, R5, R,,,, and R95 by comparing the mass arrival times of the 5%, mean, and 95% mass breakthrough fractions in the conservative and reactive transport run. The deviation of these retardation factors from the surface retardation factor (4) is a measure for the degree of nonuniform retardation of the solute plume. In addition, the ratio Deff/DL was computed for any given transport distance to estimate the relative magnitude of enhanced dispersion due to surface sorption. Figure 2.7 shows the ratio Deff/DL (upper panel) and the three 33 retardation factors (lower panel) as a function of transport distance. The ratio Deff/DL is nearly constant indicating that the extent of enhanced dispersion is independent of transport distance. The observed values of Deff/DL agree well with the independent estimate (Deff/DL=9.8, dashed line in Fig. 2.7) using the appropriate dispersion factor from Figure 2.6 (cf=14.5). Note that for transport distances on the order of the entrance length (L*crit=0.25m) , the observed ratios of Deff/DL have more uncertainty due to the poorer fit of the advectiondispersion equation to the simulated breakthrough curves, explaining the fluctuations in the data points. Although the magnitude of enhanced dispersion is nearly constant, the degree of nonuniform retardation decreases significantly with increasing transport distance (lower panel, Fig. 2.7). For example, for L=0.01m the 5% breakthrough fraction experiences only =60% of the mean retardation (Ra=5) whereas the 95% breakthrough fraction is retarded by a factor of =1.32 more than predicted by the retardation factor (Fig. 2.7). At L=0.1m the 5% breakthrough fraction experiences =90% of the mean retardation (=110% for the 95% breakthrough fraction). An additional tenfold increase in L yields essentially uniform retardation, and the three retardation factors are within 5% of the expected retardation predicted by the surface retardation factor (4). Note that the retardation of the mean arrival time is always predicted accurately by (4), independent of transport distance. 34 The decrease in the degree of nonuniform retardation is a result of the steepening of the solute breakthrough curve at greater transport distance in non-dimensional time. At small transport distances dispersive transport is a significant component of transport and enhanced dispersion has a strong influence on the retardation of the individual reference breakthrough fractions. In contrast, at long transport distances the non-dimensional BTCs are very steep, that is, the advective transport component dominates over the dispersive transport component. Although the magnitude of enhanced dispersion due to sorption (expressed by Deff) remains constant it becomes an insignificant component relative to advective transport. Hence, for long transport distances the sorbing solute is essentially transported as an "advective front". The reduction in the travel velocity of this advective front relative to the mean fluid velocity is then described by the surface retardation factor (4) . The minimum transport distance required to develop a steep solute front of uniform retardation is a function of the degree of enhanced dispersion due to surface sorption expressed by the ratio Deff/DL. In general, the required transport distances are greatest in the transport regime of Taylor dispersion which exhibit a maximum in Deff/DL (Fig. 2.5) . Similarly, a greater sorption strength (i.e. more retardation) requires a longer transport distance to approximate uniform retardation as a result of more enhanced dispersion (Fig. 35 2.6). The required transport distances vary greatly depending on which transport problem is considered. For example, for the transport problem shown in Figure 2.7, a transport distance on the order of one meter is required to approximate uniform retardation. In contrast, a strongly-sorbing solute in the transport regime of Taylor dispersion approximates uniform retardation at a transport distance on the order of several tens of meters. 2.6 Implications for Solute Transport in Rough-Walled Fractures In our analysis of the influence of surface sorption on solute dispersion we assumed that the fracture can be approximated by two parallel plates. Most natural fracture planes, however, have rough and irregular surfaces. It is uncertain to what extent our results, namely the enhanced dispersion due to surface sorption, apply to rough and/or partially closed fractures. In principle, surface roughness or partial closure would tend to disturb the uniform flow field that exists in a smooth-walled fracture of constant aperture. These perturbations of the flow field can be expected to result in additional transverse mixing of the solute by advection, and thus, they may reduce the amount of enhanced dispersion due to surface sorption. In addition, significant 36 variations in fracture aperture may not allow the development of a parabolic velocity profile in the transverse direction. The absence of this velocity profile would further reduce the amount of enhanced dispersion due to surface sorption. There is some indication that solute dispersion in roughwalled fractures is dominated by dispersion mechanisms other than Taylor dispersion. Using migration experiments in constant-aperture fractures machined from plexiglass, Dronfield and Silliman (1993) estimated dispersion coefficients for different degrees of fracture roughness. Fracture roughness was introduced by use of blockages within the fracture and/or addition of surface roughness at the fracture walls. For each fracture roughness dispersion was proportional to fluid velocity raised to an exponent. The value of the exponent, however, depended strongly on the fracture roughness, ranging from 2.0 for smooth parallel plates down to approximately 1.3 for rough plates (Dronfield and Silliman, 1993). Similarly, several laboratory migration experiments using rough-walled fractures in granite have indicated that the longitudinal dispersion coefficient is approximately proportional to the mean fluid velocity (Neretnieks et al., 1982/ Vandergraaf, pers. coram.). These experimental results suggest that rough fracture planes may behave like a porous medium rather than an open, constant aperture fracture. The reduced dependence of DL on u observed in rough fracture 37 planes suggests that (advective) mixing at junctions connecting different channels is the dominant dispersion mechanism rather than transverse diffusion across the velocity profile (Bear, 1972). In this situation, enhanced dispersion due to limited and biased sorption at the fracture walls may be insignificant. The enhanced dispersion due to surface sorption may be significant, however, in fractures where flow and transport is dominated by a few open channels. The channeling of flow and transport can be pronounced in solute transport through rough fractures in low-permeability rock (e.g., Neretnieks et al., 1982; Vandergraaf et al., 1995). Enhanced longitudinal spreading due to surface sorption is more likely to occur in flow channels, in particular if flow channeling occurs in physical channels of large aperture and/or high flow rates. It has been shown that both of these conditions tend to result in incomplete transverse mixing of the solute, and hence, tend to enhance dispersion. 2.7 Implications for Modeling Solute Transport at the Network Scale Recently published discrete network models account for sorption onto the fracture walls by introducing the surface retardation factor (4) for all fracture segments making up the network. Our results indicate that the surface retardation 38 factor provides an accurate description of the reduction in the mean tracer velocity relative to the mean fluid velocity. At the same time, our results suggest that surface sorption may result in enhanced dispersion. When using an advectiondispersion equation for transport in individual fractures (e.g., Sudicky and McLaren (1992)) the additional dispersion can be accounted for by using an effective dispersion coefficient for a sorbing solute. Our definition of the effective dispersion coefficient (equation (21) ) could be used assuming the fractures making up the network can be approximated by parallel plates. When the surface retardation factor (4) is used in network models which consider only advective transport (Dverstorp et al., 1992/ Wels and Smith, 1994), it is implicitly assumed that: (i) dispersion is negligible relative to advection and (ii) retardation is uniform in all fracture segments. The presence of enhanced longitudinal spreading of a solute plume at the scale of a single fracture (as demonstrated here for the case of a parallel-plate fracture) would limit the application of an advection-based transport model. On one hand, the magnitude of dispersion in each fracture segment would be greater than for the conservative case, in particular for a strongly-sorbing solute, thus potentially compromising assumption (i) for a greater number of fractures. On the other hand, the assumption of uniform retardation could be compromised, in particular in a dense fracture network where 39 the distances between fracture intersections can be very short. We are not aware of any studies which have focused on the influence of dispersion at the single fracture scale on dispersion at the network scale for the case of a sorbing solute. Based on our results we anticipate that the influence of dispersion at the single fracture scale on dispersion at the network scale may be considerably greater in the presence of surface sorption. This analysis would be further complicated if the network geometry introduces additional nonuniform and/or anisotropic retardation at the network scale (Wels and Smith, 1994). A consideration of nonuniform retardation at the scale of a single fracture (Fig. 2.7) may have an influence on nonuniform retardation at the network scale as well. 2.8 Conclusions Solute transport in a parallel-plate fracture has been simulated using a random walk model which accounts explicitly for sorption onto the fracture walls. Numerical simulations show that the one-dimensional advection-dispersion equation accurately describes solute transport assuming that the condition of transverse homogenization is met. A sorbing solute requires a greater transport distance for transverse homogenization to be achieved than a non-sorbing solute. This 40 increase in the required entrance length is proportional to the surface retardation factor (4). A modified entrance length has been proposed to account for this increase in the entrance length. Numerical simulations have shown that the surface retardation factor (4) gives an accurate description of the retardation of the mean tracer velocity. However, surface sorption results in enhanced longitudinal spreading of the solute plume. This additional dispersion is caused by (i) limited sorption of the solute plume to the fracture wall due to slow and incomplete transverse mixing, and (ii) biased sorption of slow-moving solute due to the parabolic velocity profile. Based on numerical simulations an effective dispersion coefficient (equation (21)) has been developed which describes this enhanced dispersion for a wide range of fluid velocities. It is shown that dispersion due to surface sorption is directly proportional to the magnitude of Taylor dispersion in a parallel-plate fracture. The increase in dispersion is negligible at very low fluid velocities and reaches a maximum in the range of fluid velocities characteristic of Taylor dispersion. It is further shown that the importance of enhanced dispersion is a function of the sorption strength of the solute, that is, an increase in retardation also results in more enhanced dispersion of the sorbing solute. At short transport distances, surface sorption results in a 41 high degree of nonuniform retardation. For longer transport distances retardation is essentially uniform since the advective component of transport dominates over the dispersive component. Nonuniform retardation may be an important consideration in the planning and interpretation of migration experiments at the laboratory scale. This effect may also be significant in the retardation of a sorbing solute at the network scale in a medium that is densely fractured. 42 2.9 References Andersson, J. and B. Dverstorp, Conditional simulations of fluid flow in three-dimensional networks of discrete fractures, Water Resour. Res. 23(10), 1876-1886, 1987. Aris, R., On the dispersion of a solute in a fluid flowing through a tube, Proc. Soc. R. Soc. London, Ser., A, 235, 6777, 1956. Bear, J., The dynamics of fluids in porous media. 7 64 pp. Elsevier, New York, 1972. Clemo, T., Dual permeability modeling of fractured media, Ph.D. thesis, Univ. of B.C., Vancouver, 1994. Dronfield, D.G. and S.E. Silliman, Velocity dependence of dispersion for transport through a single fracture of variable roughness, Water Resour. Res. 29(10): 3477-3483, 1993. Endo, H.K., J.C.S. Long, C.R. Wilson, and P.A. Witherspoon, A model for investigating mechanical transport in fracture networks, Water Resour. Res. 20(10), 1390-1400, 1984. Holtta, P., M. Hakanen, and A. Hautojarvi, Migration of radionuclides in fracture columns, In: Scientific basis for nuclear waste management XIV (T. Arbajano, Jr., L.H. Johnson, eds.) Mater. Res. Soc. P r o c , 212, 669-676, 1991. Hull, C.H., J.D. Miller, and T.M. Clemo, Laboratory and simulation studies of solute transport in fracture networks, Water Resour. Res., 23, 1505-1513, 1987. Kessler, J.H. and J.R. Hunt, Dissolved and colloidal contaminant transport in a partially clogged fracture, Water Resourc. Res., 30, 1195-1206, 1994. Kinzelbach, W. The random walk method in pollutant transport simulation, supporting paper presented in "Advances in analytical and numerical groundwater flow and quality modelling", LNE Lisbon, Portugal, June 1987. Long, J.C.S. and P.A. Witherspoon, The relationship of the degree of interconnection to permeability of fracture networks, J. Geophys. Res., 90(B4), 3087-3098, 1985. Neretnieks, I., T. Eriksen, and P. Tahtinen, Tracer movement in a single fissure in granitic rock:Some experimental results and their interpretation, Water Resour. Res., 18, 849-858, 1982. Ogata, A. and R.B. Banks, A solution of the differential 43 equation of longitudinal dispersion in porous media, Professional paper 411-A, U.S. Geol. Survey, VA, pp. A1-A7, 1961. Schwartz, F.W., L. Smith, and A.S. Crowe, A stochastic analysis of macroscopic dispersion in fractured media, Resour. Res. 19(5), 1253-1265, 1983. Water Sudicky, E.A., and R.G. McLaren, The Laplace transform Galerkin Technique for large-scale simulations of mass transport in discretely fractured porous formations, Water Resour. Res. 28(2), 499-514, 1992. Tang, D.H., E.O. Frind, and E.A. Sudicky, Contaminant transport in fractured porous media: Analytical solution for a single fracture, Water Resour. Res., 17, 555-564, 1981. Taylor, G.T., Dispersion of soluble matter in solvent flowing slowly through a tube, Proc. of the Royal Soc. of London, Series A, 219: 186-203, 1953. Taylor, G.T., Conditions under which dispersion of a solute in a stream of solvent can be used to measure molecular diffusion, Proc. of the Royal Soc. of London, Series A, 225: 473-477, 1954. Valocchi, A.J. and H.A.M. Quinodoz, Application of the random walk method to simulate the transport of kinetically adsorbing solutes, in Groundwater Contamination, IAHS Publ. No. 185: p. 35-42, 1989. Vandergraaf, T.T., D.M. Grondin, P. Vilks, and D.J. Drew, Radionuclide migration studies in the laboratory, Can. Nuclear Soc. Conf. Proceedings; 2nd Int'1 Conference on radioactive waste management, p. 142-150, 1988. Vandergraaf, T.T., C.-K. Park, and D.J. Drew, Migration of conservative and poorly-sorbing tracers in granite fractures, Radiochimica Acta, 1995 (in press). Wels, C. and L. Smith, Retardation of sorbing solutes in fractured media, Water Resour. Res., 30(9), 2547-2563, 1994. Wels, C , Chapter Three: The influence of specific surface area on solute transport in fractures - An experimental analysis, In: Transport of sorbing solutes in fractured media - A numerical and experimental analysis of dispersion and retardation, Ph.D. thesis, Univ. of B.C., Vancouver, 1995. 44 a) Infinite parallel plate 4 ^ ^ > b) 2D Model domain c) Sorption zone n Sorbed particle Diffusive step Advective step Figure 2.1. Conceptualization of reactive solute transport in a parallel-plate fracture showing (a) a pair of infinite parallel plates, (b) the two-dimensional model domain for the transport code and (c) a close-up of the "sorption zone" in proximity of the fracture wall. 45 mixed —— | Taylor —| dispersion dispersion diffusion-controlled dispersion CD 2 Pi 10 1 0 0 0 =- -P CD O CD 100 - PH CD 10 I •+J O cti 1 10 10 «-6 10 10 10 10 10 mean fluid velocity (m/s) Figure 2.2. The fracture Peclet number as a function of mean fluid velocity for b=4 00)im, Dm=2xl0~9 m2s-1 and for different transport distances L. The fracture Peclet number is not defined for fluid velocities too high to establish transverse homogenization for a given transport distance (dashed part of curves). The vertical dotted lines delineate the three transport regimes (see text). 46 1 T * >°° 0.8 ^ /P o u " 0.6 • u = 5 x l 0 ms ^ u = l x l u ms O u=oxlO ms 0.4 0.2 0 0.6 0.8 1 1.2 1.4 N o n - d i m e n s i o n a l time (ut/L) Figure 2.3. Simulated breakthrough curves of a conservative solute for mean fluid velocities representing (i) diffusioncontrolled dispersion (circles), (ii) mixed dispersion (triangles), and (iii) Taylor dispersion (solid dots). All other transport parameters are constant (L=0.2m, b=400|im, Dm=2xl0~9 m2s_1) . The lines show the advection-dispersion equation using the respective fracture Peclet numbers Pef=20 (dotted line), 3450 (solid line), and 105 (dashed line), respectively. 47 1 i 0.8 o o o ppb^gisno^ i r D. 0" .-"O ^IO -,<cr L=0.2m 0.6 0.4 A R a =1.0 • Ra = 2.0 o Ra=51 0.2 o°y * i 0 i i i J I I L 0.5 1 1.5 2 2.5 N o n - d i m e n s i o n a l time (ut/LR a ) Figure 2.4. Simulated breakthrough curves of (i) a conservative solute (Ra=l), (ii) a weakly-sorbing solute (Ra=2), and (iii) a strongly-sorbing solute (Ra=51). The conservative transport parameters are the same as those used for Taylor dispersion in Figure 2.3 (u=5xl0~3ms_1, L=0.2m, b=400jim, Dra=2xl0~9 it^s-1) . The best-fit Peclet number for the weakly- and strongly-sorbing solute are Pe=14.0 (dashed line) and Pe=3.7 (dotted line), respectively. The fracture Peclet number for the conservative BTC is Pef=105 (solid line) . 48 diffusion—controlled dispersion mixed —— | Taylor —| dispersion dispersion 8 6 Q 0 10 8 10 ? 10 6 10 5 10 4 10 3 10" 2 mean fluid velocity (m/s) Figure 2.5. The ratio Deff/DL as a function of the mean fluid velocity for a weakly-sorbing solute (Ra=2) . The solid dots represent results from individual numerical migration experiments. The solid line shows the expected behaviour (equation (21)). The vertical dotted lines delineate the three transport regimes. 49 25 1 1—i—i—i i i i 1 1—i—i—r-TT J I I •4-H o u o o S o w 5M r—1 CD PH 20 15 10 5 00 Q J 0' I I I I I I 1 10 1 I I I I 100 Retardation factor (Rf) Figure 2.6. The dispersion factor cf determined numerically for a range of retardation factors. The value of cf is independent of the physical parameters (u,b,DJ used in the simulations. 50 11 — i — i — i i i 111 I I I Mill i i i i 111 R a =5.0 10.5 Q m 10 CD 9.5 9 0.01 10 i i i i i 111 j 1—I—I 1 I II i i i i i 11 1 r & Rm A R5 " 7.5 K L 1 0.1 1 J f -H—B- f^f i i ii ii I 2.5. 0 0.01 J I II II 1 0.1 Figure 2.7. The ?Jb?o 1 § e P f HJ t (u$jir : %SS?) ^ ^ r e t a r d a t i o n factors Rref for the 5%, mean , and 95% breakthrough fractions (lower panel) as a function of transport distance for a sorbing solute (Ra=5). The required entrance length for the sorbing solute is LRa>0.25m. 51 Chapter Three: THE INFLUENCE OF SPECIFIC SURFACE AREA ON SOLUTE TRANSPORT IN FRACTURES AN EXPERIMENTAL ANALYSIS 3.1 Introduction A number of experimental studies have been performed to elucidate the processes that control reactive solute transport in low-permeability fractured media. Neretnieks and co-workers carried out migration experiments using the sorbing radionuclides strontium and cesium in natural, rough-walled fractures in granite (Neretnieks et al., 1982/ Moreno et al., 1985). Their results suggested that reactive solute interacts with the granite in two ways. First, it sorbs directly to the fracture walls (surface sorption) causing retardation of the advective movement in the fracture. Second, the solute diffuses into the matrix where it sorbs at inter- and intragranular surfaces (matrix sorption), resulting in a loss of solute from the fracture and strong tailing of the solute breakthrough curve. Their interpretation of the reactive breakthrough curves was, however, complicated by the rough walls of the fracture and the dispersion that occurred due to channeling. Neretnieks et al. (1982) modeled the observed 52 dispersion assuming the solute traveled in a set of independent, parallel channels with different channel apertures. Neretnieks et al. pointed out that the (advective) movement of the sorbing tracer displayed a stronger dependence on the channel aperture (tracer velocity proportional to aperture cubed) than did a non-reacting tracer (fluid velocity proportional to aperture squared) . Hence, a detailed description of the range of apertures in the fracture plane (or channel apertures in the channel model) appears crucial when describing reactive transport in a rough-walled fracture. In a follow-up experiment, Moreno et al. (1985) demonstrated that either an advection-dispersion model or their channel model could be fit to the experimental breakthrough curves. However, the choice of a particular dispersion model influenced the parameter estimates describing surface sorption and matrix diffusion/sorption. The interpretation of reactive transport experiments in fractures is further complicated when the solute exhibits nonideal sorption behaviour or when chemical heterogeneity is present. Both conditions are known to influence solute dispersion in porous media (e.g. Brusseau and Rao, 198 9; Tompson, 1993) . Chemical heterogeneity is usually defined as spatial heterogeneity in sorption sites and/or sorption strength of the sorbing material. Non-ideal sorption behaviour includes (i) non-linear sorption, (ii) chemical nonequilibrium, and (iii) hysteresis in sorption. Hysteresis 53 between adsorption and desorption is frequently observed, in particular for strongly-bound ions, with desorption taking more time than adsorption to attain equilibrium (Davis and Kent, 1990). In many cases hysteresis results in stronger sorption, that is, a greater partitioning of solute mass between the solid phase and solution, during desorption compared to during adsorption. The influence of chemical heterogeneity and non-ideal sorption on solute dispersion in fractures is best examined in a constant aperture fracture. With this experimental set-up, Taylor dispersion is the only physical dispersion mechanism. Taylor dispersion of conservative and ideally-sorbing solutes in a constant-aperture fracture has been described using analytical techniques (Kessler and Hunt, 1994) and numerical simulations (Wels, 1995) . Some migration experiments of sorbing solutes have been performed in constant aperture fractures (Vandergraaf et al., 1988/ Fujikawa et al., 1993). Vandergraaf et al. (1988) observed for a system with continuous tracer injection, that the dispersion of cesium in sawed fractures in granite was considerably greater than that of a conservative solute. Vandergraaf et al. noted that deviations from the linearity and reversibility of the sorption process could have resulted in front sharpening or tailing of the elution profile. A quantitative analysis of the breakthrough curves was not presented. Fujikawa et al. (1993) performed similar migration 54 experiments using Dirac pulse injections of cesium in sawed and chiseled fractures in granite. They observed a strong dependence of the cesium breakthrough curves, in particular peak height and tailing, on the fluid velocity in the fracture. A transport model that accounted for surface sorption and matrix diffusion/sorption was fit to the experimental data. Fujikawa et al. found that transport parameters characterizing the sorption process had to be varied significantly to obtain a good fit for the different experimental breakthrough curves of cesium. The authors proposed that microscopic heterogeneities in the diffusion/sorption potential of the rock matrix could be responsible for the observed discrepancies in parameter estimates. The authors noted that a consideration of the nonlinear sorption may have improved the model fit. The possibility of non-equilibrium sorption and/or hysteresis in sorption was not considered. The influence of non-ideal sorption behaviour on solute transport is difficult to assess in the case of rough-walled fracture planes. Neretnieks et al. (1982) suggested that the non-linear nature of the cesium isotherm (as determined from batch experiments) could have caused the discrepancies they observed between peak outflow concentrations and model predictions based on a linear sorption model. Vandergraaf et al. (1995) injected several radionuclides of varying sorption strengths into a rough-walled fracture (0.8m by 0.9m) in a 55 quarried block of granite. The movement of the tracers was inferred on the basis of breakthrough curves from individual boreholes and by scanning the sorbed activity on the fracture surfaces. Conservative breakthrough curves suggested that the solute travelled predominantly in an open channel with an aperture on the order of 1mm. The peak retardation of the weakly-sorbing strontium was consistent with surface sorption coefficients determined in static sorption experiments assuming transport had occurred in the dominant channel (Vandergraaf et al., 1995). However, the moderately-sorbing cesium appeared to have travelled in two "separate phases". While a small fraction of the injected cesium was eluted, the bulk of the cesium mass remained sorbed within the dominant channel of the fracture (Vandergraaf, 1995). This complex transport behaviour could be a result of physical heterogeneity causing channeling, chemical heterogeneity, nonideal sorption behaviour, or some combination of the above. None of these experimental studies has focused specifically on the influence of specific surface area on transport of a sorbing solute in a fracture. The specific surface area, aw, of a fracture is the ratio of the fracture surface area to volume of mobile water. To illustrate the influence of aw consider the transport equation for a sorbing solute in a single fracture. Assuming fast, linear, and reversible sorption on the fracture surface and within the rock matrix, the differential equation for solute transport in the fracture 56 is ( e . g . Moreno e t a l . , dCf dt 1985) L d2Cf dx2 dCf dx ^ e dcj dz z=0 (1) where De is the effective diffusion coefficient of the matrix as defined by Neretnieks (1980) (De= 8pT Dm) , and Ra is the surface retardation factor, defined as Ra = 1 + awKa (2) All other parameters are defined in the notation. Assuming only diffusive mass transfer in the rock matrix the transport equation for solute movement within the porous matrix is ^ dt = ^ Kv ^ dz2 (3) where Kv is the so-called volumetric distribution coefficient. It is related to the sorption capacity of the rock KDps by (Neretnieks, 1980) Kv = ep+KDps 57 (4) where KD is the bulk distribution coefficient. The specific surface area, aw, influences solute transport in two ways. First, aw enters into the definition of the surface retardation factor (equation (2)). A larger aw implies a greater surface area for surface sorption to take place relative to the amount of tracer in solution. As a result, a greater amount of tracer is sorbed, the surface retardation factor increases, and the solute plume in the fracture advances at a slower rate. Second, the specific surface area also influences the fourth term in equation (1) which describes the rate of exchange of the solute with the matrix by molecular diffusion. This term implies that a larger specific surface area provides more access of the solute to diffuse from the open fracture into the rock matrix relative to the amount of tracer in solution. Note that these relationships also hold if the fracture volume varies while the surface area is held constant. In short, a sorbing solute travelling in a fracture with a high specific surface area (i.e. a narrow aperture) should experience more retardation and more tailing than the same solute travelling in a wider fracture, assuming rock properties are the same. Note that the influence of specific surface area on solute transport depends on the sorption strength of the solute (characterized by Ka and K D ) . The stronger the solute sorbs 58 onto the fracture walls and within the matrix, the larger will be the influence of aw on solute transport in the fracture. For a non-reactive solute (Ka=KD=0) the influence of aw is negligible unless the porosity is relatively large resulting in significant storage of the solute within the matrix (see equation (4)). The influence of specific surface area on reactive solute transport at the scale of a single fracture may have important implications for reactive solute transport at the network scale. Wels and Smith (1994) performed numerical migration experiments in fracture networks assuming advective transport of a sorbing solute is retarded in each fracture segment according to equation (2). It was shown that the coupling of retardation and fracture aperture (=2/aw) at the scale of a single fracture produces nonuniform and anisotropic retardation of the solute plume at the network scale (Wels and Smith, 1994) . Moreno and Neretnieks (1993) studied the influence of aw on reactive solute transport in fractured media, considering both surface sorption and matrix diffusion/sorption. They demonstrated that the magnitude of the specific surface area has a strong impact on the retardation of the solute plume. In their study, Moreno and Neretnieks refer to aw as the "flow-wetted surface" of a fracture. We prefer the use of "specific surface area" indicating that aw is the ratio of the fracture surface area to the fracture volume with units m2/m3. 59 In this chapter we study the influence of specific surface area on the transport of the sorbing solutes strontium and uranium in smooth-walled, constant-aperture fractures. For this purpose migration experiments were performed in machined fractures which differed only in fracture aperture. The simple geometry of the fracture allowed an interpretation of experimental breakthrough curves by means of the transport equations (l)-(4). The surface and bulk distribution coefficients Ka and KD were estimated independently using static sorption experiments. The experimental results presented here demonstrate that the influence of specific surface area on solute retardation may be significantly greater than predicted by the surface retardation factor (2) if the solute shows hysteresis in the sorption process. 3.2 3.2.1 Materials and Methods Migration Experiments The migration experiments were performed in rectangular, constant-aperture fractures in steel and granite. Our analysis focuses on the transport in granite. Steel was chosen for comparison. The use of steel minimizes the effects of matrix diffusion and spatial heterogeneity in sorption potential on solute transport. The fractures were constructed by mounting 60 two slabs of granite (or steel) together (total dimensions -10.0x2.5x2.5 cm3) . The fractures in steel were obtained by milling a 2.3cm wide slot along the length of one slab to the desired depth (=aperture). Prior to assemblage, the surfaces of both steel slabs were oxidized by placing the pieces for 72 hours in an oven at 500°C. This procedure provided a smooth, fine-grained and uniform surface coating of iron oxides (magnetite with some hematite) which served as a sorbing substrate for the reactive tracer. The granite slabs were cut from quarried rock blocks of the Lac du Bonnet batholith, Manitoba, Canada (Vandergraaf et al., 1982). The fractures were made by placing narrow teflon spacers of desired thickness (=aperture) between two granite slabs. The surfaces of the slabs facing the fracture were polished to minimize the roughness of the fracture walls. The outer surfaces of the mounted slabs were sealed using silicone rubber to avoid evaporation from the rock matrix. Two different experimental set-ups were used in the migration experiments. For a continuous or step input of tracer, an inflow reservoir was used and solutions tagged with radionuclides were passed through the fracture using a gravity feed system (Fig. 3.1a). A teflon reservoir (=5x5x5cm3) was mounted directly on top of the joined slabs allowing the solution to enter the open fracture directly through a narrow slit (-2.3x0.lcm2 ) in the base of the reservoir. The outflow from the fracture was funneled via a triangular slit cut into 61 a teflon plenum into narrow teflon tubing (I.D.=0.025cm). A tubing length of 0.9 m provided sufficient resistance to obtain the desired flow rates. The background solution in the inflow reservoir could be manually replaced by tracer solution or vice versa in a short time period (<2 min). For a pulse injection, tracer was delivered using a High Pressure Liquid Chromatography (HPLC) pump. In this set-up both ends of the joined slabs were mounted directly to teflon plenums which connected to inflow and outflow tubing (Fig. 3.1b). The tracer solution was injected as a "slug" into the inflow tubing just upstream of the fracture using a syringe (Fig. 3.1b). In both set-ups the outflow from the core was collected in a fraction collector, after which the samples were analyzed for the radioactive tracers (U-233/3H20 or Sr90/3H2O) using liquid scintillation counting. A total of five migration experiments was performed using the sorbing radionuclide U-233 in the fractures made of steel and the radionuclide Sr-90 in the fractures made in granite (Table 3.1). To study the effect of specific surface area, migration experiments were run in a narrow and wide fracture for each material (Table 3.1). The experiments also differed with respect to the injection scheme used (Table 3.1). In experiments 1 and 2 the tracer solution was injected continuously for the entire duration of the experiment. In experiment 3 the tracer solution was injected as a short pulse (duration of the injection, At=0.35hr). In experiments 4 and 5 62 the tracer solution was injected as a long step-input (At=25 and 46hr, respectively). In all cases, the non-reactive tracer tritium (3H20) was injected simultaneously with the sorbing tracer to estimate physical transport parameters. The chemical composition of the tracer input solutions is summarized in Table 3.2. A consistent loss of tritium on the order of 10-13% was observed in all migration experiments. This loss is believed to be a result of evaporation during sampling. These evaporative losses seem plausible given the small sample volumes involved (1-4 drops=0.05-0.2 ml) and the experimental conditions (flow rates on the order of 20 drops/h, low humidity, fume hood). The tritium breakthrough curves were adjusted to account for this loss. The small sample volumes did not allow a monitoring of pH in the outflow. Some change in solution pH during the migration experiments can thus not be ruled out. Fluctuations in solution pH would complicate the interpretation of the uranium migration experiments. Uranium is known to form several hydroxide species for solution pH£ 5.5 which may differ in their sorption behaviour. Possible pH fluctuations during migration experiments are not believed to influence the sorption characteristics of strontium since it does not form any hydroxide complexes for solution pH^ 10. After termination of the migration experiments the cores were dried and disassembled. Autoradiographs of the fracture 63 surfaces were then taken in order to estimate the spatial distribution of the radioactive tracer remaining sorbed on the fracture surfaces. 3.2.2 Static Sorption Experiments Batch sorption measurements were carried out for Sr-90 with crushed granite in order to determine the bulk sorption capacity of the granite matrix. The rock sample consisted of crushed Lac du Bonnet granite (e.g. Vandergraaf et al., 1982) which was wet-sieved through a 35-80 mesh (particle diameters ranging from 180-435 Jim) . Samples of one gram of dried rock were placed in polyethylene bottles and contacted with 10ml of tracer solution containing Sr-90. Five initial tracer concentrations were used ranging from 50-2000 Bq/ml. All samples were shaken manually every 15 minutes and allowed to equilibrate for a total of eight hours. Small aliquots of contacting solution were taken after 1, 4 and 8 hours and counted for strontium. The strontium concentrations remained constant in control runs containing no rock samples indicating that sorption to the walls of the bottle was negligible; hence all loss of tracer from solution can be attributed to sorption. Coupon sorption experiments were performed to obtain independent estimates of the sorption capacity of the 64 macrosurfaces of the polished granite and the oxidized steel. The experimental procedure was similar to that described by Vandergraaf and Abry (1982). Small coupons were cut from core material, saturated, and all sides but the one with the fracture surface were sealed with silicone rubber. The fracture surfaces were contacted with tracer solution by pipetting solution (2-3cm3) directly onto the surface (4-6cm2) . A narrow rim of silicone rubber around the edges of the fracture surface prevented the solution from spilling over. In a sorption step, tracer solution with a chemical composition representative of that used in the migration experiments was added to the surface and allowed to equilibrate. The equilibration time for sorption ranged from 1-24 hours. In a subsequent desorption step, the tracer solution was replaced by a tracer-free solution. The equilibration time for desorption was the same as previously allowed for sorption of a given coupon. The coupons (and contacting solution) were gently shaken throughout the experiment. In a first set of experiments the contacting solutions were sampled only once, that is after the end of the specified equilibration time for a given coupon. The change in tracer concentration in the solution during contact was used to calculate a surface distribution coefficient for the given equilibration time representing either sorption (Kas) or desorption (Kad) . These Ka values represent bulk estimates of solute partitioning which may include effects of matrix 65 diffusion/sorption, in particular for large equilibration times. In a second set of experiments (granite coupons only) small aliquots of the contacting solution were taken repeatedly to measure the change in tracer concentration over time. A diffusion-sorption model was fit to these time trends of tracer concentration (see below). The model distinguishes between (instantaneous) surface sorption and (time dependent) matrix dif fusion/sorption. The values of Kas (or Kad) obtained from a model fit identify the tracer partitioning due to "instantaneous" surface sorption (desorption) only. In all coupon sorption experiments, evaporation rates were measured to allow a correction of tracer concentration due to evaporation of solution. 3.2.3 Interpretative Models The tracer breakthrough curves were interpreted using an analytical solution to the transport problem introduced earlier (equations (1-4)). Note that the tracer concentrations are measured in the column outflow, and as such they represent a "flux-concentration" (e.g. Kreft and Zuber, 1978). Assuming zero initial concentrations in the fracture and matrix, a constant concentration at the inlet of the fracture during tracer injection, and a matrix and fracture of very large 66 extension, given the flux concentration, ( e . g . Moreno e t a l . , Cf(x,t) in the f r a c t u r e 1985) ^f 2 ,Pe.r (PetwG/8z2) Pe x r _ , 22 PPe* e2 _ — = — e x p ( — - ) e x p ( - r 2 - - — - )erfc ){ [t-(Pet0/4r2)] 1 6 r c l o y/fr is ) dz (5) i where Pe = ux/DL (6) I = (Pet0/4t)1/2 (7) t 0 = i ? a x / u = J? a t w (8) G = a w (D e K v ) 1 / 2 (9) Note that the parameter G describes the influence of matrix diffusion/sorption on the concentration distribution in the fracture. This parameter is a function of the diffusion/sorption capacity of the matrix (DeKv) and the specific surface area of the fracture. Since the migration experiments were performed in smooth-walled, constant-aperture fractures, the specific surface area is approximated by two times the inverse of the fracture aperture (i.e. a w =2/b). According to Fujikawa and Fukui (1990), a flux boundary at 67 the inlet (and a "flux injection" during injection time) is a more appropriate boundary condition for fracture column experiments than a constant concentration boundary used above. A comparison of the analytical solutions for a flux-type boundary condition (eq. 2-7 in Fujikawa and Fukui (1990)) and for a constant concentration boundary (see equation (5) above) showed that the two solutions are equivalent in advectiondominated systems (Pe>10) (see also Maloszweski and Zuber, 1990). The transport conditions in our experiments were always advection-dominated (Pe>10). For computational convenience we prefer the use of equation (5) in the breakthrough curve analysis. Since the transport problem (1-4) is linear in Cf superposition in time applies and the flux concentration for a finite step injection is obtained by subtracting two Cf solutions shifted in time by the duration of the step injection At (Kreft and Zuber, 1978) . The second set of sorption experiments on granite coupons (see above) was interpreted using a "static" diffusionsorption model that describes the change in tracer concentration in the contacting solution due to surface sorption and matrix diffusion/sorption. The transport equation for matrix diffusion/sorption is given by equation (3). The tracer concentration in the contacting solution as a function of time is given by equation (1) by setting u and DL equal to zero. In these coupon experiments, the specific surface area aw is given as the inverse of the height h of the water column 68 contacting the coupon. In a static coupon sorption experiment, injection and detection are in the resident fluid (i.e. in the contacting solution) and the delta-type fluid injection takes the form (Fujikawa and Fukui, 1990) Co - J - (10) where V is the volume of the contacting solution, and M is the total mass (or activity) of tracer added. Note that the retardation factor Ra appears in the initial condition (10). Under the assumption of instantaneous sorption, the solute sorbs to the surface at the moment of tracer injection and therefore, from the beginning, the solute exists in both the solution phase and on the surface of the coupon. The principle of mass conservation requires that the initial tracer concentration in the contacting solution be reduced by the magnitude of the retardation factor Ra (see Fujikawa and Fukui (1990) for details) . Assuming zero initial concentration throughout the rock matrix and diffusion distances much shorter than the thickness of the granite coupons (i.e. <<2cm), the change in tracer concentration in the contacting solution is given by (Neretnieks, 1980) 69 c -f = exp (6) erfc (6 i n ) (11) where the dimensionless time 0 is given as e = ^x (12) t h2 The diffusion-sorption model is used to determine Ka values describing instantaneous surface sorption as follows. First, the analytical solution (11) is fit to the observed concentration-time trends using a fixed value of the DeKv parameter estimated from batch sorption measurements of KD and literature values of De# ps and £p. In this way, the initial concentration C0 defined in (10) is the only fitting parameter. Second, the surface distribution coefficient Ka is calculated from the best-fit C0 and the (known) mass M of tracer initially introduced according to K^M^-U (13, The above equation was obtained by re-writing the initial condition (10) in terms of Ka for t = 0. Note again that this 70 approach assumes that sorption to the surface of the coupon occurs instantaneously. The calculated Ka values express the loss of tracer from solution due to instantaneous sorption only. Tracer loss due to matrix diffusion/sorption is accounted for explicitly in the model by the diffusion/sorption parameter DeKv. 3.3 Results and Discussion 3.3.1 Influence of Fracture Aperture on Retardation of Uranium Migration experiments 1 and 2 indicated a strong retardation of U-233, presumably due to strong sorption of uranium onto the oxidized steel surfaces of the fracture. For example, in migration experiment 1 the uranium concentration in the outflow was less than one percent of C0 even after the flushing of 54 fracture volumes. The high retardation of uranium made a detailed analysis of the breakthrough curve infeasible. Instead, the retardation of uranium was estimated from the spatial distribution of the tracer sorbed onto the fracture surfaces. Figure 3.2 shows the spatial distribution of the uranium plume in the wide (top) and narrow (bottom) fractures as captured by the autoradiograph of the fracture surfaces. As 71 expected, the uranium plume had advanced further in the wider fracture than in the narrow fracture (Fig. 3.2). Based on the autoradiograph we estimated an average plume advance of 3.5cm and 1.6cm in the wide and narrow fracture, respectively. Accounting for the duration of the experiment (54 and 45 fracture volumes in experiment 1 and 2, respectively) the bulk retardation of the uranium plume was approximately Ra=154 in the wide fracture (800jim) and Ra=2 82 in the narrow fracture (400|im) . In other words, a twofold decrease in fracture aperture resulted in an almost twofold increase in retardation. The observed difference in retardation agrees fairly well with the theory. According to equation (2) the retardation of the solute front is inversely proportional to fracture aperture for a strongly-sorbing (Ra>>10) solute such as uranium. Static sorption experiments of U-233 on oxidized steel coupons suggested an average value for the surface distribution coefficient of Ka~l. 0 ±0.3cm (n=8 samples). Using this "static" estimate of Ka we obtain a retardation factor of Ra=26 ±8.5 for the wide fracture and Ra=51 ±15 for the narrow fracture. These values are approximately six times lower than the retardation factors observed in the migration experiments. In our opinion differences in the pH and the redox conditions of the static solution and the transport solution could explain these discrepancies. Although the starting compositions of transport and static solutions were very 72 similar, they may have changed during the experiment. For example, we observed the formation of amorphous iron-oxide on the coupon surfaces after several hours of contact during static sorption. In contrast, the oxidized steel surfaces of the fractures remained (macroscopically) unaltered, even several days after starting the experiment. A better control on the experimental conditions (e.g. use of an anoxic chamber and buffer solutions) would probably result in a better agreement of Ra values estimated from static and migration experiments. An attempt was made to obtain the concentration profile of the uranium sorbed onto the fracture surfaces. The oxide coating was scraped off in narrow increments of 0.25cm and leached in 6N HC1 to redissolve the sorbed uranium. However, only ~ 8 % of the total mass of uranium injected was recovered using this procedure. Apparently, the majority of the uranium mass had diffused through the removable layer of the oxide coating and was sorbed directly onto the steel surface. The concentration profile of the "removable uranium" is shown in Figure 3.3. The influence of fracture aperture is also recognizable in the retardation of this uranium fraction. Using spatial moments, the centre of mass appeared to have moved twice as far in the wide fracture compared to the narrow fracture. This observation, in conjunction with the overall spatial distribution of the uranium plume (Fig. 3.2), provides qualitative support for the concept that retardation of a 73 strongly-sorbing solute is inversely proportional to fracture aperture. 3.3.2 Static Sorption Experiments of Strontium on Granite Batch sorption measurements of Sr-90 on crushed granite suggested a linear Freundlich isotherm (exponent=0.99; R2=0.98) for the concentration range of interest (Fig. 3.4). The influence of equilibration time (1-8 hrs) on the distribution of tracer between rock surface and solution was negligible. The least-squares fit of the Freundlich isotherm to all data points suggested a bulk distribution coefficient of KD=29.5 ml/g (Fig. 3.4). The linearity of the strontium sorption isotherm has been demonstrated in other batch sorption experiments (e.g. Skagius et al., 1982) . However, the bulk distribution coefficient determined in our experiment is higher than what has been observed in other batch sorption experiments of strontium on crushed granite samples. Torstenfelt et al. (1982) reported KD values ranging from 5-20 ml/g for strontium on granite samples from Sweden. Vandergraaf et al. (1995) determined a KD value of only 1.3 ±0.1 ml/g for strontium on similar granite material as used in our experiments. The low ionic strength of the background solution used in this study (Table 3.2) is likely responsible for the relatively high bulk distribution coefficient observed here. A 74 lower ionic strength of the background solution provides fewer cations to compete with strontium for the existing sorption sites. As a result, a greater amount of strontium sorbs to the solid phase and the value of KD increases. The surface distribution coefficients for strontium, estimated from the coupon sorption experiments on granite, are summarized in Table 3.3. Recall that a "bulk estimate" of Ka is based on the total partitioning of the solute between the solid phase and solution (including possible sorption within the matrix) . The "model estimates" of Ka shown in Table 3.3 are discussed in conjunction with the diffusion-sorption model (see below). The variability in "bulk estimates" of Ka was relatively small in comparison to differences in Ka as a result of equilibration time or direction (sorption versus desorption) (Table 3.3). For example, the average surface distribution coefficient determined in a sorption step for short contact times (1-2 hrs) was Kas=0.07 ±0.03cm compared to Kas=0.19 ±0.02cm for six hours contact time. Longer contact times also resulted in larger differences between the surface distribution coefficients determined during sorption (Kas) and desorption (K/) (Table 3.3). The influence of equilibration time on sorption was studied in more detail in an additional experiment. Three granite coupons were exposed to tracer solution and the change in tracer concentration monitored over time (Fig. 3.5). The mass (activity) of strontium initially introduced to the solution 75 is shown for comparison (as concentration M/V, Fig. 3.5). Figure 3.5 shows that most of the sorption took place during the first hour (or possibly minutes) of contact, after which sorption increased at a fairly slow rate. The sorption time trends were very similar for all three coupons (Fig. 3.5). Torstenfelt et al. (1982) observed similar time trends of strontium sorption on granite macrosurfaces. They suggested that the first rapid phase corresponds to sorption on externally accessible surfaces, while the second phase might reflect a diffusion-controlled mass transport into microfissures and internal exchange sites. In order to separate the effects of fast surface sorption and slow matrix diffusion/sorption the observed concentrationversus-time trends are modeled using the analytical solution (11) of a diffusion-sorption model. Recall that in this model the parameter Ka describes instantaneous surface sorption and the parameter DeKv describes matrix diffusion/sorption. For simplicity, we use an independent estimate of DeKv and fit the model by varying only C0. Assuming an effective diffusivity De=l. OxlO-12 m2/s as determined experimentally by Skagius et al. (1982) for crushed granite particles, and further using the volumetric distribution coefficient Kv=78.2 based on KD=29.5 ml/g determined from the batch sorption isotherm (Fig. 3.4), the dif fusion/sorption parameter becomes DeKv=7 . 8xl0-11 m 2 /s. This value of DeKv is an order of magnitude higher than DeKv values determined experimentally by Moreno et al. (1985). This 76 difference is expected and can be attributed to a higher sorption strength (higher KD) observed in our system (see above). The best fits of the diffusion-sorption model to the data are shown in Figure 3.5 as solid lines. The best-fit value of C0 is indicated for reference. The diffusion-sorption model fits the experimental data well, suggesting that matrix diffusion/sorption could explain the slow decline in tracer concentration in solution over time (Fig. 3.5). The "model estimate" of Ka is calculated from the best-fit value of C0 and the (known) mass M of tracer initially introduced (see equation (13)). For the three granite coupons shown in Figure 3.5 we obtain model estimates of Kas (describing instantaneous surface sorption) ranging from 0.09cm to 0.13cm (Table 3.3). Note that the estimates of Kas would be significantly higher (2-3 times) if matrix diffusion/sorption was not considered explicitly and instead Kas values were estimated from the total partitioning after 4 hours (so-called "bulk estimates"). Clearly, the use of the diffusion-sorption model provides a more detailed description of surface sorption, in particular in those cases where matrix diffusion/sorption is significant. The desorption trends as a function of time for the three granite coupons are shown in Figure 3.6. The mass (activity) of strontium present in solution prior to initial desorption was negligible. An initial, fast desorption was followed by a relatively slow rate of further release of tracer into 77 solution with time. Unfortunately, the diffusion-sorption model could not be applied here due to uncertainty in the initial and boundary conditions. As a first approximation, a straight line was fit to the experimental data to obtain a rough estimate of the amount of strontium desorbed "instantaneously" (Fig. 3.6). This estimate was then used to calculate surface distribution coefficients Kad representing instantaneous desorption for any given coupon. According to model calculations a significant fraction of the tracer (-60%) had moved into the matrix during the sorption step. Assuming that this tracer fraction did not participate in the initial desorption the calculated "model estimates" of Kad ranged from 0.16cm to 0.48cm (Table 3.3). Significantly higher surface distribution coefficients were obtained (Kad=0.9-1.3cm) when the influence of matrix diffusion/sorption was not considered, that is by calculating the total partitioning of strontium between surface/matrix and solution after four hours desorption. A comparison of the various estimates of Kas and Kad (Table 3.3) shows that the Kas values are systematically lower than Kad values suggesting hysteresis in the sorption process. In the case of bulk estimates, some of the discrepancy between Kas and Kad is caused by matrix diffusion/sorption which releases a sorbed tracer only slowly from the matrix ("apparent hysteresis"). However, the effects of matrix diffusion/sorption are filtered out in the model estimates of 78 Kas and Kad. The fact that the model estimate of Kad is significantly greater than the estimate of Kas strongly suggests that strontium exhibited hysteresis during sorption onto granite macrosurfaces under the experimental conditions used here. Hysteresis between adsorption and desorption of ions is frequently observed, with desorption taking more time than adsorption to attain equilibrium (Davis and Kent, 1990). Batch experiments by Barney (1984) demonstrated that the magnitude of this hysteresis depends not only on the solute chemistry but also on the experimental conditions (e.g., temperature, redox conditions, and ionic strength). The low ionic strength of the background solution used in our experiments (Table 3.2) probably favoured slow desorption and/or partiallyirreversible sorption as few ions could compete with strontium for the sorption sites. Vandergraaf and Abry (1982) observed an even higher degree of hysteresis for Sr-90 sorption/desorption on granite during coupon sorption experiments using synthetic groundwater and longer equilibration times (14-28 days). In their experiments surface distribution coefficients for sorption (Kas=0.09-0.17cm) and desorption (Kad=l.8-2.6cm) differed by one order of magnitude (Vandergraaf and Abry, 1982). Unfortunately, the reported values of Ka are bulk estimates only. It is not clear how much of this discrepancy is due to apparent hysteresis caused by matrix diffusion/sorption. 79 In summary, the interpretation of coupon sorption experiments with a diffusion-sorption model has shown that strontium exhibited significant hysteresis between adsorption onto and desorption from granite macrosurfaces. The influence of this hysteresis on the transport of strontium in fractures in granite is discussed in later sections. 3.3.3 Analysis of Tritium Breakthrough Curves The breakthrough curves (BTCs) of the non-sorbing tracer tritium were used to determine the hydraulic properties (Pe, tw) during the migration experiments in granite. The relatively short residence times (tw~l hr) did not allow a determination of the G parameter (equation (9)) which takes into account the interaction with the rock matrix. For a nonsorbing solute this parameter takes the form G=aw (Deep) 1/2. Given the range of values one can expect for the effective diffusivity of tritiated water in granite (De=0 . 7-1. 8xl0~13 m 2 /s; Moreno et al., 1985) and porosity in Lac du Bonnet granites (ep=0.02-0.04), the influence of G on the breakthrough curve is negligible. We assumed De=1.0xl0~13 m2/s and £p=0.03 in order to determine the remaining parameters (Pe, tw) . Figure 3.7 shows the experimental BTCs of tritium observed in experiment 3 (top) and 5 (bottom), respectively. The nearsymmetrical shape of the BTCs suggests that a Fickian-type 80 dispersion model is applicable. The analytical solution fits the data fairly well (Fig. 3.7). Significant outliers in the data are believed to be a result of experimental error such as variability in drop size during sampling. The hydraulic parameters Pe and tw are summarized in Table 3.4. From these values the dispersion coefficient DL and the total volume Vt are calculated. The total volume Vt includes "dead volume" from tubing and inflow/outflow plenums and is calculated from the measured flow rate and the residence time tw. The values of Vt agree to within 10% with independent estimates of Vt based on the calculated volumes of fracture and tubing. The relatively high Peclet numbers (>100) indicate that advective transport dominated over dispersive transport in all experiments (Table 3.4). A comparison of the DL values indicates that the amount of hydrodynamic dispersion was similar in all experiments (Table 3.4). The observed solute dispersion is approximately an order of magnitude greater than expected for Taylor dispersion across the fracture aperture alone (e.g., Kessler and Hunt, 1994/ Wels, 1995). Dye experiments in acrylic replicates of the fractures suggested that a velocity profile across the width of the fracture (see Fig. 3.1a) contributed significantly to hydrodynamic dispersion in our experiments. The latter velocity profile is believed to be an entrance effect introduced by the inflow reservoir. Apparently, the experimental set-up used for the pulse 81 injection resulted in more dispersion than the set-up used for the continuous injection (Table 3.4). It is likely that some mixing occurred between the injected tracer slug and the replaced background solution before the tracer slug entered the fracture. This "smearing" of the tracer slug could explain the greater dispersion in the tritium elution profile in experiment 3. However, the experimental set-up (e.g. tracer injection, tubing) probably contributed much less to the observed dispersion of the reactive solute strontium. In the following section we will show that the interaction of strontium with the fracture surface and the granite matrix dominated the dispersion of strontium. 3.3.4 Influence of Non-Ideal Sorption on Dispersion of Strontium In migration experiment 3 we focus on the influence of nonideal sorption on the transport of strontium in a fracture. Figure 3.8 shows the observed breakthrough of strontium (symbols) in the outflow in response to a short pulse injection. The upper panel (Fig. 3.8a) shows the instantaneous breakthrough curve as a function of dimensionless effluent volume. A dimensionless effluent volume of one corresponds to an effluent volume equal to one fracture volume Vt (including dead volume). The observed mass recovery is shown, in the lower 82 panel (Fig. 3.8b). The mass recovery curve illustrates how much of the injected tracer mass is lost from the fracture due to matrix diffusion/sorption and/or irreversible sorption (e.g. Maloszewski and Zuber, 1990). The experimental breakthrough curve of strontium shows a steep rise with an early peak and a gradual decline with significant tailing (Fig. 3.8a). This asymmetry in the breakthrough curve was not observed for tritium (Fig. 3.7a). Clearly, strontium exhibits more dispersion than the non-reactive tracer. Note that this "enhanced dispersion" results in lower retardation of the peak (Rpeak=6.4) relative to retardation of the centre-of-mass (R50=10.0). Using the definition of the surface retardation factor (2), we obtain surface distribution coefficients ranging from Ka=0.21cm (based on peak retardation) to Ka=0.35cm (based on centre-of-mass retardation). These estimates of the sorption strength fall within the range of Ka values determined in coupon sorption experiments (Table 3.3). We now use the transport model (1-4) to estimate the degree of enhanced dispersion of strontium. Note that the model accounts explicitly for the dispersion and tailing caused by matrix diffusion/sorption. The enhanced dispersion due to sorption is modeled implicitly by increasing the value of the dispersion coefficient beyond the value of DL determined for the non-reactive tracer tritium. By using this "effective" dispersion coefficient (Deff) we assume that enhanced dispersion due to sorption is "Fickian". We will compare three 83 scenarios with (i) hydrodynamic dispersion only (Deff=DL=0 .7 6cm2/s; Pe=100) , (ii) enhanced dispersion (Deff=7. 6cm2/s; Pe=10) , and (iii) strong enhanced dispersion (Deff=38cm2/s; Pe=2) . We adopt the same diffusion/sorption capacity DeKv=7 . 8xl0~n m2/s as used earlier for modeling the coupon sorption experiments (see above). The value of Ka which describes the sorption potential to the fracture surface is held constant in all three simulations. We use the average value of the model estimates (Table 3.3) for Kas (instantaneous sorption) and Kad (instantaneous desorption), that is Ka=0.25cm. We emphasize that the model is used here for illustrative purposes rather than for fitting of transport parameters. The results of the three transport simulations (lines) are shown in Figure 3.8. An inspection of Figure 3.8a shows that the transport simulation with Pe=100 (based on hydrodynamic dispersion alone) does not fit the observed breakthrough of strontium. A Peclet number Pe=100 produces a narrower, more symmetric peak, and significantly higher peak concentrations than are actually observed (Fig. 3.8a). Only the tail of the breakthrough curve is reproduced fairly well by this model run (Fig. 3.8a). In other words, the process of matrix diffusion/sorption may explain the significant tailing. However, it can not explain the observed broadening of the peak and the reduction in peak height. Clearly, additional dispersion has occurred as a result of chemical interaction 84 with the fracture surfaces. The good fit of the transport model with the observed strontium data for Pe=10 suggests that the dispersion of strontium was approximately an order of magnitude greater than hydrodynamic dispersion. An inspection of the simulated recovery curves (Fig. 3.8b) shows that the "final" mass recovery is virtually independent of the assumed degree of solute dispersion in the fracture. Instead, diffusion into and sorption within the rock matrix controls the tailing behaviour of the recovery curves and the final mass recovery. Using our independent estimate of the matrix dif fusion/sorption parameter (DeKv=7 . 8xl0-11 m2/s) the transport model predicts a tracer mass recovery of =85% after 40 dimensionless effluent volumes (Fig. 3.8b). These predictions agree fairly well with the observed strontium recovery of =92% suggesting that our estimate of DeKv (based on KD from the batch isotherm) is a good first approximation for the diffusion/sorption parameter. The following three processes may have contributed to the additional dispersion observed for strontium: (i) chemical heterogeneity, (ii) hysteresis in sorption, and (iii) limited transverse mixing across the fracture aperture. The influence of chemical heterogeneity on dispersion of sorbing solutes is well documented for the case of a porous medium (e.g. Brusseau and Rao, 1989; Tompson, 1993). The same principles would apply to a fracture, namely an increase in solute dispersion with an increase in spatial variability in sorption sites and/or 85 sorption strength. Autoradiographs of the granite fracture surfaces, taken after completion of experiment 3, revealed a very heterogeneous ("spotty") spatial distribution of the strontium remaining sorbed on the surface. This chemical heterogeneity is at the scale of the mineral grains (l-2mm). Similar spatial distributions of sorbed activity were found in experiments 4 and 5. These observations are consistent with results from several sorption studies which have demonstrated that the various minerals making up a granite vary in their general sorptive capacity as well as their specific affinity for sorption of strontium (e.g. Torstenfelt et al., 1982/ Ticknor et al., 1986). Hence, chemical heterogeneity could have contributed to the enhanced dispersion of strontium. Hysteresis of sorption may also cause enhanced dispersion due to differences in retardation of the front and tail of the tracer pulse. Retardation of the front of the tracer pulse is controlled by sorption to the fracture walls whereas the retardation of the tail of the tracer pulse is controlled by desorption from the fracture walls. It follows that the tail of the pulse is retarded more (higher Kad) than the front of the pulse (lower Kas) causing an increase in plume spreading. The significant hysteresis of strontium observed during "static" sorption to granite macrosurfaces (Kad > Kas, Table 3.3) probably contributed to the enhanced dispersion of strontium in the migration experiment. The third process that could cause enhanced dispersion, 86 limited transverse mixing, relates to the limited access of the tracer to the sorbing surface during advection-dominated transport in a constant aperture fracture (Wels, 1995). The parabolic velocity profile which develops across the fracture aperture results in fast advective transport in the central region of the fracture whereas solute in the proximity of the fracture walls experiences little or no advective movement. At the same time, sorption of the solute occurs only at the fracture walls. The degree to which the solute plume participates in surface sorption depends on the degree of transverse mixing across the aperture. The degree of mixing in a constant-aperture fracture is characterized by the transverse Peclet number Pet=ub/Dm. For Pet<2 the diffusion distance b is sufficiently small to approximate "instantaneous transverse mixing" (Wels, 1995). The above condition implies that the entire solute plume participates in surface sorption. For Pet>2 instantaneous transverse mixing is not approximated and transverse mixing is "limited". In the case of a fast, reversible sorption, limited transverse mixing results in enhanced dispersion of the solute plume (Wels, 1995). In experiment 3, the transverse Peclet number was Pet=14.8, suggesting that transverse mixing was indeed limited. According to numerical results by Wels (1995), the effective dispersion coefficient for strontium (assuming Ra=7.4) is approximately eight times greater than the longitudinal dispersion coefficient of a conservative solute under the 87 transport conditions in experiment 3. However, these numerical results assume a parallel-plate fracture (no dispersion across the fracture width) and ideal sorption. Both assumptions are not fully satisfied in our experiment which makes an accurate estimation of the magnitude of enhanced dispersion due to limited transverse mixing difficult. We can only speculate on the relative importance of chemical heterogeneity, hysteresis, and limited transverse mixing in causing dispersion of strontium during experiment 3. The coupon sorption experiments suggest that chemical heterogeneity and hysteresis are interrelated. The range of "instantaneous" Ka values determined for three different coupons (i.e. chemical heterogeneity) was much greater during the desorption step (Fig. 3.6) than during the sorption step (Fig. 3.5). More in-depth studies on these complex, non-ideal sorption processes are needed in order to quantify their respective influence on solute dispersion in fractures. In summary, the sorbing solute strontium exhibited approximately an order of magnitude more dispersion than predicted from the breakthrough of the non-reactive tracer tritium. This enhanced dispersion is believed to be a combination of chemical heterogeneity, hysteresis in the sorption process, and limited transverse mixing. 88 3.3.5 Influence of Fracture Aperture on Retardation of Strontium The influence of fracture aperture on strontium transport is best illustrated by direct comparison of the breakthrough curves of strontium in experiments 4 and 5. Figure 3.9 shows the outflow concentrations of strontium as a function of dimensionless effluent volume. In experiment 4, two short, but significant deviations from the general trend were observed as a result of sudden changes in strontium concentrations in the inflow reservoir. The spike of C/C0 observed at about six effluent volumes (Fig. 3.9) was in response to an (accidental) increase of 10% in the strontium concentration in the inflow reservoir. The short but significant drop in C/C0 at about 18 effluent volumes (Fig. 3.9) was in response to switching from strontium containing input solution to tracer-free background solution (end of step injection). These short-term effects are not considered further in our analysis. Two general observations can be made regarding the influence of fracture aperture on strontium transport. First, the retardation of strontium is significantly greater in the narrow fracture than in the wide fracture (Fig. 3.9). For example, the first significant breakthrough of strontium (C/C0>0.025) was observed after 1.5 fracture volumes in the wide fracture. In the narrow fracture, however, approximately 8.7 fracture volumes were required to elute the same 89 breakthrough fraction (Fig. 3.9). Second, the strontium breakthrough is much more dispersed and tracer concentrations are lower in the fracture with narrow aperture compared to that with a wide aperture. For example, the peak effluent concentration in the wide-aperture fracture was approximately Cpeak/Co=0.55 compared to only Cpeak/Co=0.30 in the narrowaperture fracture. These observations agree qualitatively with the theoretical predictions from the transport model (1-4). According to the model a larger specific surface area (i.e. a smaller aperture) results in (i) increased retardation and (ii) increased loss of tracer from the fracture to the matrix and tailing in the effluent. In the following, we compare the experimental data to simulated breakthrough curves using the transport model (14 ) . The purpose of this comparison is to test whether the differences in strontium breakthrough in experiments 4 and 5 can be explained by the difference in fracture aperture. In theory, a single set of parameter values for Pe, Ka and DeKv should describe both experimental breakthrough curves. An extensive comparison of transport simulations with the experimental data for a wide range of Pe, Ka and DeKv values yields two general conclusions. First, the transport model fit the breakthrough curves in both experiments very poorly with Pe=250 which is the estimate of Pe from the tritium BTC (Table 3.4). The Peclet number had to be decreased by an order of magnitude to account for enhanced dispersion due to non-ideal 90 effects of sorption. Second, no single parameter set (Pe, Ka, and DeKv) could be found that would describe the experimental data in both experiments equally well. In particular, Ka had to be varied greatly to obtain acceptable fits of the model simulations with the data. To illustrate this point we will show three model "type curves" for experiments 4 and 5, assuming the following values of the surface distribution coefficient: (i) Ka=0.1cm, (ii) Ka=0.5cm, and (iii) Ka=1.0cm. For simplicity, we set Pe=10 which assumes that the magnitude of enhanced dispersion in experiments 4 and 5 is equal to that estimated for experiment 3. In these model calculations we also assume DeKv-7. 8xlO-11 m2/s as used earlier for modeling the sorption time trends (Fig. 3.5) . Figures 3.10 and 3.11 show the simulated type curves (lines) for experiments 4 and 5 together with the experimental data (symbols). Instantaneous breakthrough curves are shown in the upper panel and mass recovery curves are shown in the lower panel. In experiment 4, the model type curve with Ka=0.lcm describes the retardation in the arrival time of the peak of the plume fairly well (Fig. 3.10a). This value of Ka is within the range of model estimates of Ka determined from coupon sorption experiments (Table 3.3). In contrast, the same type curve (Ka=0.1cm) significantly underestimates the retardation of strontium observed in experiment 5 (Fig. 3.11a) . The Ka value had to be increased by one order of magnitude (Ka=l.0cm) to model the observed retardation in this experiment. 91 These model results indicate that strontium sorbed more strongly to the fracture walls during experiment 5 than during experiment 4. Additional results not shown here indicated that strontium retardation was much more sensitive to Ka than to DeKv, suggesting that differences in surface sorption rather than in matrix diffusion/sorption caused the large differences in retardation between the two experiments. Note that the difference between the Ka values required to match the strontium breakthrough curves is much greater than the variability in the Ka values determined in coupon sorption experiments (Table 3.3). Hence, we do not believe that natural variability in sorption strength between the fracture surfaces caused the different retardation response in experiments 4 and 5. In our opinion, the difference in sorption strength of strontium is caused by the difference in specific surface area of the fracture in experiment 4 (b=7 80jim) and experiment 5 (b=450jim) . This contention is supported by the fact that the simulation results for experiment 4 are consistent with those for experiment 3 which was also conducted in a wide fracture (b=780|im) . In both experiments, the estimates of sorption strength agree fairly well (Ka=0.25cm in experiment 3 versus Ka=0.1cm in experiment 4) despite the differences in tracer injection (Table 3.1). We emphasize that our results suggest that the distribution of tracer concentration between surface and solution (i.e. Ka) is a function of specific surface area 92 aw. This is not to be confused with the commonly assumed influence of aw on the retardation factor Ra ( = l+awKa) which follows from the partitioning of tracer mass between surface and solution (expressed by awKa) . A positive correlation of Ka with aw results in a much stronger (non-linear) influence of specific surface area on retardation (as observed in our experiments) than is commonly assumed in the retardation model (equation 2 ) . Two different factors could have contributed to the apparent increase in sorption strength with an increase in specific surface area aw. First, it is conceivable that the number of sorption sites on the fracture walls is not infinite (as assumed in the Freundlich isotherm model). In the case of a limited number of sorption sites, a higher specific surface area in a fracture implies more sorption sites and thus a greater potential for surface sorption resulting in a higher distribution coefficient. In theory, this hypothesis could be tested by determining Ka values in coupon sorption experiments for a series of surface area-to-volume ratios (=SA/V=aw) . Vandergraaf and Drew (1991) used such an approach to study the influence of specific surface area on the sorption of cesium to granite coupons. They found that a doubling of the surface area-to-volume ratio (from SA/V=0. 8cm"1 to 1.6cm"1) resulted in an increase in Ka estimates of up to 100%. Our coupon sorption experiments did not suggest a systematic dependence of the strontium sorption on specific surface area for ratios ranging 93 from SA/V«l-4cm-1. Note, however, that the SA/V ratios used in the coupon experiments were approximately an order of magnitude lower than the specific surface areas of the fractures used in the migration experiments (e.g. aw=25.6cm~1 for b=7 80]nm) . The second factor which could have influenced the apparent sorption strength relates to the limited access of the tracer to the sorbing surface during advection-dominated transport in a constant-aperture fracture (see previous section). Recall that a transverse Peclet number greater than two indicates that mixing across the aperture is limited due to high advection relative to transverse diffusion. The transverse Peclet numbers for experiment 4 (b=7 80jim) and experiment 5 (b=450jim) were Pet=10.0 and Pet=5.7, respectively. In other words, transverse mixing and access to the sorbing fracture walls was limited to a greater extent in the wide fracture (greater diffusion distance) relative to the narrow fracture. In the following we hypothesize that the difference in the degree of transverse mixing, in conjunction with hysteresis in sorption, could explain the different sorption strength of strontium in experiments 4 and 5. First, consider the transport of strontium in the wide fracture. In experiment 4, the final tracer recovery (Fig. 3.10b) was significantly lower (-60%) than predicted by the type curve for Ka=0.1cm which describes the retardation of the eluted strontium. Apparently, the extent of desorption of 94 strontium was so small (i.e. Kad >> Kas=0. lcm) that a significant fraction of the injected strontium sorbed "irreversibly". We hypothesize that limited transverse mixing promotes irreversible sorption by limiting the exchange of sorbed tracer with dissolved tracer moving by advection in the central region of the fracture. This hypothesis is consistent with the apparently weak sorption strength (Ka=0.lcm) of the eluted strontium. Much of the eluted strontium probably traveled in the central region of the fracture with little opportunity to sorb to the fracture walls. Note that matrix diffusion/sorption is an unlikely sink for the strontium not recovered. Model calculations suggest that the DeKv value would have to be raised by more than an order of magnitude (>8xl0~10 m2/s) to reduce the predicted final tracer recovery to near 60%. In our opinion, partially irreversible sorption as a result of limited transverse mixing is a more likely cause for the low tracer recovery than excessive matrix diffusion/sorption. Next, consider the case of strontium transport in the narrow fracture which allows greater diffusive mixing across the aperture. In experiment 5, the trend in tracer recovery (Fig. 3.11b) is simulated fairly well with the type curve for Ka=1.0cm. Here, the observed tracer recovery is consistent with the (albeit very high) retardation of the eluted strontium. According to our hypothesis, the increased transverse mixing in experiment 5 resulted in participation of 95 most of the injected strontium in surface sorption. It follows that desorption affected most of the eluted strontium which explains the very high Ka value required to model the retardation of the breakthrough curve (Ka=l.Ocm, see Fig. 3.11a). In short, the greater the degree of transverse mixing, the greater will be the effect of hysteresis (high Kad) on retardation of the eluted strontium. Note that hysteresis had a much stronger influence on strontium transport in the case of a finite step injection on the order of 1-2 days (experiments 4 and 5) compared to a short pulse injection on the order of less than one hour (experiment 3 ) . For example, the retardation of the solute plume in experiment 3 could be modeled using the average of the instantaneous estimates (Table 3.3) for sorption and desorption (i.e. Ka= (Kas+Kad)/2= 0.25cm) . The modeling of the strontium retardation in experiment 5, however, indicated a much stronger "average" distribution coefficient (Ka=l.Ocm). Assuming a decrease in the extent of desorption is responsible, we obtain an approximate distribution coefficient for desorption of Kad=l. 9cm (=2Ka-Kas) . In experiment 4, such an estimation of Kad from the "average retardation" of strontium is not meaningful since a large fraction of strontium sorbed irreversibly. It is clear, however, that the extent of desorption was much lower (i.e. higher Kad) than estimated for experiment 3. We believe that the different magnitude of hysteresis is a 96 result of the difference in sorption time. A longer duration of tracer injection results in a longer residence time of the tracer sorbed onto the surface. Results of coupon sorption experiments not shown here suggested that an increase in equilibration time during a sorption step decreased the extent of desorption during a subsequent desorption step. Longer sorption times may allow the formation of surface complexes (chemi-sorption) and/or the incorporation of the tracer into the mineral lattice. These new chemical bonds can be expected to have greater activation energies for dissociation (e.g. Davis and Kent, 1990). The possible influence of sorption time on hysteresis and its effects on transport of the sorbing solute is a subject worth further study. In conclusion, retardation of strontium was approximately 20 times greater during transport in a narrow fracture (b=450jim) compared to that in a wide fracture (b=780p.m) . This increase in retardation is much greater than predicted by the commonly used definition of the surface retardation factor (equation (2)). A combination of hysteresis of sorption and limited transverse mixing could explain the observed strong influence of fracture aperture (specific surface area) on the retardation and mass recovery of strontium. 97 3.3.6 Implications for Reactive Transport in Fractured Media Most authors have adopted the concept of retardation in their analysis of transport of sorbing solutes in fractured media (e.g. Neretnieks et al., 1982; Moreno et al., 1985; Dverstorp et al., 1992; Fujikawa et al., 1993; Moreno and Neretnieks, 1993; Wels and Smith, 1994). The retardation model assumes that retardation of the solute in the fracture varies in proportion to the sorption strength (Ka) and the specific surface area (aw=2/b) (see equation (2)). This model implies that retardation increases with a decrease in fracture aperture whereas the sorption strength of a solute is independent of fracture aperture (or aw) . The assumed (linear) dependence of retardation on fracture aperture has important implications for transport in fractured media. Numerical studies have shown that this dependence would result in nonuniform and anisotropic retardation at the network scale (Wels and Smith, 19 94). By analogy, the same concept may apply to reactive transport in rough-walled fractures where retardation in a given channel is inversely proportional to channel aperture (Neretnieks et al., 1982). Our experimental results suggest that the influence of fracture aperture on retardation may be greater than assumed by the retardation equation (2) for those sorbing solutes that exhibit hysteresis in surface sorption. Our findings suggest that the sorption strength (Ka) may increase with decreasing 98 fracture aperture (increasing aw) thus magnifying the influence of specific surface area on retardation. The observed dependence of Ka on fracture aperture would significantly increase the complexity of solute transport either in rough-walled fractures or in networks of fractures. In particular, we anticipate a greater degree of nonuniformity and anisotropy in retardation at the plume scale. It has been suggested that hysteresis in sorption caused the significant increase in strontium retardation in the narrow fracture. However, other non-ideal effects in sorption such as non-linear and/or non-equilibrium sorption may also cause a non-linear dependency of retardation on fracture aperture. Given its potential to significantly influence retardation at the plume scale, the relationship between retardation and fracture aperture should be examined for a range of solutes and transport conditions. It should be noted that migration experiments in natural, rough-walled fractures may not be ideal -for this purpose due to the difficulties in characterizing the aperture distribution and the resulting dispersion and/or channeling. In our opinion, migration experiments in fractures of simple geometry are best suited for studying the influence of fracture aperture on retardation. 99 3.4 Conclusions Laboratory migration experiments using sorbing radionuclides (U-233 and Sr-90) in constant-aperture fractures (oxidized steel and granite, respectively) were carried out to study the influence of specific surface area on reactive solute transport. The radionuclides were injected in fractures of different apertures while maintaining similar fluid velocities (fluid residence times ~1 hour). Uranium was found to sorb strongly to the oxidized surfaces of the fractures in steel. Autoradiographs of the uranium plume sorbed to the fracture surfaces and spatial moments of "removable uranium" extracted from the oxide coating suggest that the uranium had advanced approximately twice as far in the wide fracture (800]um) compared to that in the narrow fracture (400jim) . These observations provided qualitative support for the concept that retardation of a strongly-sorbing solute is inversely proportional to fracture aperture. Strontium exhibited non-ideal sorption behaviour on granite surfaces under the experimental conditions studied here. Static sorption experiments on granite coupons suggested a hysteresis in the sorption process, showing higher surface distribution coefficients for desorption than for sorption (i.e. Kad > Kas) . A diffusion-sorption model which accounts explicitly for diffusion into and sorption within the rock matrix was used to estimate Ka values describing 100 "instantaneous" surface sorption. These model estimates ranged from Kas=0 . 09-0 .13cm for instantaneous sorption to Kad=0.160.48cm for instantaneous desorption. In the migration experiments, strontium experienced significantly more dispersion than the non-reactive tracer tritium. This enhanced dispersion is believed to be a result of chemical heterogeneity, hysteresis in sorption, and limited transverse mixing across the fracture aperture. Model simulations suggested that the dispersion of strontium was approximately an order of magnitude greater than that due to hydrodynamic dispersion. The observed retardation of the peak (Rpeak=6.4) and the centre-of-mass (R50=10.0) of the strontium breakthrough curve was consistent with model simulations using an independent estimate of the surface distribution coefficient from static sorption experiments (Ka=0.25cm). The influence of fracture aperture on retardation was much greater than predicted by the commonly used definition of the surface retardation factor (equation (2)). Strontium retardation was approximately an order of magnitude greater in a narrow fracture (450jim, Ra~45) than in a wide fracture (780pm, Ra~3.5) using step injections of one and two days, respectively. We hypothesize that hysteresis in sorption, in conjunction with limited transverse mixing across the aperture, caused the apparent increase in sorption strength with a decrease in fracture aperture. The experimental design had a significant influence on the 101 mass recovery of strontium. For example, 92% of the injected strontium was recovered in response to a short pulse injection (At=0.35 hr) compared to only 60% in response to a finite step injection (At=24 h r ) . These results suggest that the magnitude of hysteresis increases with the residence time of the sorbing solute in the fracture. 102 3.5 Notation aw specific surface area of fracture, cm 2 /cm 3 b fracture aperture, jam Cf concentration in the mobile fluid in the fracture, Bq/ml Cp concentration in the stagnant fluid in the matrix, Bq/ml Cr concentration in solution contacting rock coupons, Bq/ml C0 initial concentration, Bq/ml Dm molecular diffusion coefficient in water, cm 2 /h De effective diffusivity of rock matrix, £ptDm, cm 2 /h DL longitudinal dispersion coefficient, cm 2 /h G parameter defined in (9) h height of water column in static sorption experiments, cm KD bulk distribution coefficient, ml/g Kv volumetric distribution coefficient, cm 3 /cm 3 Ka surface distribution coefficient, cm Kas surf, distr. coeff. determined in sorption step, cm Kad surf, distr. coeff. determined in desorption step, cm 1 parameter defined in (7) M Activity of tracer injected into resident solution, Bq Pe Peclet number, ux/D L Q volumetric discharge, ml/h Ra surface retardation factor t time, h t0 tracer residence time, h tw water residence time, h 103 u fluid velocity in fracture, cm/s V volume of solution contacting rock coupons, ml Vf volume of fracture, ml Vt total volume of fracture and tubing, ml x distance along the fracture, m z distance into rock matrix, m ps rock density, g/cm3 £p porosity of rock matrix X tortuosity of the rock matrix 0 parameter defined in At time duration of tracer injection, h 3.6 (12) Acknowledgements We thank Dr. T.T. Vandergraaf and the staff of the Applied Geochemistry Branch of the Whiteshell Nuclear Establishment for their generous support during this study. Doug Drew provided important technical support and Kenneth V. Ticknor prepared the autoradiographs. Dana L. Joseph from the Analytical Science Branch performed most of the radiometric analyses. This research has also been supported by a grant to L. Smith from the Natural Science and Engineering Research Council of Canada (NSERC). 104 3.7 References Barney, G.S., Radionuclide sorption and desorption reactions with interbed materials from the Columbia River basalt formation, Proceedings of the 185th American Chemical Society Meeting, Div. of Nuclear Chem. and Tech., p. 3-23, 1984. Brusseau, M.L. and P.S.C. Rao, Sorption non-ideality during organic contaminant transport in porous media, Crit. Rev. in Env. Control 19(1), 33-99, 1989. Davis, J.A. and D.B. Kent, Surface complexation modeling in aqueous geochemistry, in "Mineral-water interface Geochemistry" Hochella and White (eds.), Reviews in Mineralogy, 23, Mineral Soc. of America, 1990. Dverstorp, B., J. Andersson, and W. Nordqvist, Discrete fracture network interpretation of field tracer migration in sparsely fractured rock, Water Resour. Res. 28, 2327-2343, 1992. Fujikawa, Y., M. Fukui, D.J. Drew, and T.T. Vandergraaf, Analysis of the migration of instantaneously injected cesium in artificial fractures of Lac du Bonnet granite, Manitoba, Canada, J. Contam. Hydrol., 14, 207-232, 1993. Fujikawa, Y. and M. Fukui, Adsorptive solute transport in fractured rock: analytical solutions for delta-type source conditions, J. Contam. Hydrol., 6, 85-102, 1990. Kessler, J.H. and J.R. Hunt, Dissolved and colloidal contaminant transport in a partially clogged fracture, Water Resourc. Res., 30, 1195-1206, 1994. Kreft, A. and A. Zuber, On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions, Chem. Eng. Sci., 31, 1471-1480, 1978. Maloszweski, P. and A. Zuber, Mathematical modeling of tracer behaviour in short-term experiments in fissured rocks, Water Resourc. Res., 1517-1528, 1990. Moreno, L., I. Neretnieks, and T. Eriksen, Analysis of some laboratory tracer runs in natural fissures, Water Resour. Res., 21, 951-958, 1985. Moreno, L. and I. Neretnieks, Flow and nuclide transport in fractured media: The importance of the flow-wetted surface for radionuclide migration, J. Contam. Hydrol., 13, 49-71, 1993. Neretnieks, I., Diffusion in the rock matrix: an important factor in radionuclide retardation? J. Geoph. Res., 85, 4379105 4397, 1980. Neretnieks, I., T. Eriksen, and P. Tahtinen, Tracer movement in a single fissure in granitic rock: Some experimental results and their interpretation, Water Resour. Res., 18, 849858, 1982. Skagius, K., G. Svedberg, and I. Neretnieks, A study of strontium and cesium sorption on granite, Nucl. Techn., 59, 302-313, 1982. Tang, D.H., E.O. Frind, and E.A. Sudicky, Contaminant transport in fractured porous media: Analytical solution for a single fracture, Water Resour. Res., 17, 555-564, 1981. Ticknor, K.V., T.T. Vandergraaf, and D.C. Kamineni, Radionuclide sorption on mineral and rock thin sections, Part II: Sorption on granitic rock, A.E.C.L. Technical Report TR 385, 1986. Tompson, A.F.B., Numerical simulation of chemical migration in physically and chemically heterogeneous porous media, Water Resour. Res., 29(11), 3709-3726, 1993. Torstenfelt, B., K. Andersson, and B. Allard, Sorption of strontium and cesium on rocks and minerals, Chem. Geol., 36, 123-137, 1982. Vandergraaf, T.T., Abry, D.R.M., and C.E. Davis, The use of autoradiography in determining the distribution of radionuclides sorbed on thin sections of plutonic rocks from the Canadian Shield, Chemical Geology 36, 139-148, 1982. Vandergraaf, T.T. and D.R.M. Abry, Radionuclide sorption on drill core material from the Canadian Shield, Nucl. Technol., 57, 399-412, 1982. Vandergraaf, T.T. and D.J. Drew, Laboratory studies of radionuclide transport in fractured crystalline rock, Proceedings of the 3rd International Symposium on Advanced Nuclear Energy Research - Global Environment and Nuclear Energy, Mito City, Japan, 1991, pp. 385-393, 1991. Vandergraaf, T.T., D.M. Grondin, and D.J. Drew, Laboratory radionuclide migration experiments at a scale of one meter, Mater. Res. Soc. Symp. Proc., 112, 159-168, 1988. Vandergraaf, T.T., C.-K. Park, and D.J. Drew, Migration of conservative and poorly-sorbing tracers in granite fractures, Radiochimica Acta, 1995 (in press). Vandergraaf, T.T., Radionuclide migration experiments under 106 laboratory conditions, Geophysical Research Letters, (in press) Wels, C. and L. Smith, Retardation of sorbing solutes in fractured media, Water Resourc. Res., 30(9) , 2547-2563, 1994. Wels, C , Chapter Two: The influence of surface sorption on solute dispersion in parallel-plate fractures, In: Transport of sorbing solutes in fractured media: A numerical and experimental analysis of dispersion and retardation, Ph.D. thesis, Univ. of B.C., Vancouver, 1995. 107 Table 3.1. 3xp. # Experimental conditions in migration experiments. Tracers Material Injection At(h) b(jim) V£(ml) Q(ml/h) 1 233 U/ 3 H steel continuous 50 800 1.88 1.71 2 233 U/ 3 H steel continuous 38 400 0.94 0.78 pulse 0.35 780 1.44 1.47 3 90 Sr/ 3 H granite 4 90 Sr/ 3 H granite finite step 25 780 1.44 1.06 5 90 Sr/ 3 H granite finite step 46 450 0.97 0.67 108 Table 3.2. 3xp. 1 2 # Chemical composition of tracer solutions. Tracer Input C0 (Bq/ml) Background PH 233 u ( 3 H 2 0) 14 . 0 (188) 10" 3 M NaCl 6 .1 233 u ( 3 H 2 0) 15.5 (94) 10" 3 M NaCl 6 .1 3 90 Sr ( 3 H 2 0) 430 (916) 10" 5 M HN0 3 5.5 4 90 Sr ( 3 H 2 0) 99 . 3 (617) 10" 6 M HNO3 5.7 5 90 Sr ( 3 H 2 0) 99 . 3 (617) 10" 6 M HNO3 5.7 109 Table 3.3. Summary of surface distribution coefficients determined in coupon sorption experiments of strontium on polished granite. Type of Estimate Contact time (h) Number of Samples Kas Kad cm Bulka 1-2 n=4 0.07±0.03 0.11+0.03 Bulka 6 n=3 0.19+0.02 0.32+0.01 "instantaneous" n=3 0.09-0.13 0.16-0.48 Modelb a average and st. dev. from "bulk measurements"; no correction made for matrix diffusion/sorption b range of Ka values for three coupons describing instantaneous sorption (Fig. 3.5) and instantaneous desorption (Fig. 3.6) (see text) 110 Table 3.4. Summary of hydraulic parameters determined for migration experiments in fractures of granite using the transport model (1-4) based on tritium. Exp. # Fluid1 residence time tw, h 3 1.3 1.62 100 0 .76 4 1.4 1.47 250 0.29 5 1.6 1 .07 250 0.25 Peclet1 Number Pe Total Volume Vt, ml Dispersion2 Coefficient DL, cm2/h 1 based on visual fit; uncertainty estimated to be +10% 2 calculated assuming x=10cm 111 Figure 3.1a. Experimental set-up for continuous injection. 112 Outflow plenum I Granite sample | jlnflow plem Tracer injection hose I Figure 3.1b. Experimental set-up for pulse injection. 113 (a) Experiment 1 (b=800 |am) • continuous injection for 54 fracture volumes (Vf) • av. plume advance -3.5 cm • approximate Ra= 154 (b) Experiment 2 (b=400 |am) • continuous injection for 45 fracture volumes (Vf) • av. plume advance -1.6 cm • approximate Ra= 282 F i g u r e 3 . 2 . A u t o r a d i o g r a p h s of t h e uranium plume on t h e f r a c t u r e s u r f a c e i n experiment 1 (top) and experiment 2 ( b o t t o m ) . The dark a r e a s i n d i c a t e t h e p r e s e n c e of s o r b e d a c t i v i t y of uranium. 114 20 E x p e r i m e n t 1 (b = 800/xm) 15 10 CD la IflTnnnm4rn 0 0 i 6 8 10 10 E x p e r i m e n t 2 (b=400/zm) O CO CO c\2 I 0 .N^RN 0 2 4 6 8 Distance f r o m s o u r c e (cm) 10 Figure 3.3. The concentration profile of removable uranium on the fracture surfaces in experiments 1 (top) and experiment 2 (bottom). 115 10 I I I I I I A 1 hr D 4 hrs -A 8 h r s tao u* a 10 — i — i — i i i 11 i i — i—i—i [Sr] s o r b =29.8 [Sr] s o l i i i it 0.99 = x* £ 1000 CO 100 J 1 IA I I I I I 10 100 Sr in solution (Bq/ml) I I I l lI 1000 Figure 3.4. Experimental results of batch sorption of strontium on crushed granite particles (l80-435um) for 1, 4 and 8 hrs equilibration time. A least-squares fit to all data (line) suggests a linear Freundlich isotherm (R2=0.98). 116 Sa o •H +-> o (72 CO Sorption time (hrs) Figure 3.5. Sorption of strontium to three granite coupons as a function of time. The best fit of the diffusion-sorption model (solid line) assuming 0^=7 . 8x10'1X m2/s was used to determine C0 and Kas (see text) . 117 PQ o •i-H •+-> (72 0 1 2 3 4 Desorption time (hrs) Figure 3.6. Desorption of strontium from three granite coupons as a function of time. The initial concentration of strontium in solution was negligible (stars). A straight line was fit to the data to estimate C0 and a first approximation of K,d (see text) . 118 0.8 -| 1 r "i 1 r n 0.6 i r Pe=100 0.4 o \ u 0.2 OB—'—l 0 1 0.8 0.6 0.4 0.2 0 0 1 2 3 Effluent v o l u m e (ml) Figure 3.7. Experimental breakthrough of tritium in response to a pulse injection (experiment 3, top) and a finite step injection (experiment 4, bottom). The lines represent the best (visual) fit of the transport model (1-4) . 119 1 1 1 ,1 Experiment 3i 1 1 1 j 1 1 II 1 i 1 i i i 1 _ o 0.03 o 1 :"\ \ X I—I \ \/jf\ 0.02 _;' - — —- Pe=2 — Pe=10 - - Pe=100 ; \ "(a) Jmf*\ \ 0.01 — o CO ^g 0 i 'i i 0 CD > o o K i i i i 1 1 1 1 1 i i I "i 20 10 1 0.8 1 1 1 1 i 1 1 1 I r—• 30 i L(b) - i i _ 40 i i • I • • i _ _ — - -*^™*^ 0.6 0.4 ~~ r -''j? •' _/ / M! / JF1 0.2 Jr i - / Jr I o CO 0 - / JF<* I • ^ - " 0 i 1 i 10 i 1 20 i i i i 1 i 30 Dimensionless Effluent Volume i i i 40 ut/x Figure 3.8. Experimental breakthrough of strontium (dots) and model simulations for experiment 3 assuming DeKv=7 . SxlO"11 m2/s and (i) Pe=2 (dashed line), (ii) Pe=l0 (solid line), (iii) Pe=100 (dotted line). Normalized effluent concentration is shown in the top panel and normalized mass recovery in the bottom panel. 120 0.8 o u i r i • Exp. 4 • Exp. 5 0.6 • u 0.4 • ofeffc - r • 3D D o -p w g • 0.2 • § 03 0 %•#• J 0 20 L • • • • 40 Dimensionless Effluent Volume 60 ut/x Figure 3.9. Comparison of experimental breakthrough of strontium in experiment 4 (b=7 8 0um; open dots) and experiment 5 (b=450um; solid dots). 121 Experiment 4 o o o X I—I o CO > o o a K o •+J CO 0 20 40 Dimensionless Effluent Volume 60 ut/x Figure 3.10. Experimental breakthrough of strontium (dots) and model type curves for experiment 4 (b=7 8 0um) assuming (i) Ka=0.1cm (solid line), (ii) Ka=0.5cm (dotted line), (iii) Ka=1.0cm (dashed line). The diffusion/sorption parameter is DeKv=7 . 8xl0-11 m2/s in all three simulations. 122 Experiment 5 o O o X fa a • i-H •+J o •+-> 01 S 1 1 I 1 1 1 1 1 1 1 1 I 1 1 1 \ S 0.8 L- (b) - - > O o a> K 0.6 - ^ / 0.4 - • • " •'' ' --' - — 0.2 PI o -4-5 CO 0 0 20 40 60 Dimensionless Effluent Volume 80 ut/x Figure 3.11. Experimental breakthrough of strontium (dots) and model type curves for experiment 5 (b=450um) assuming (i) Ka=0.1cm (solid line), (ii) Ka=0.5cm (dotted line), (iii) Ka=1.0cm (dashed line). The diffusion/sorption parameter is DeKv=7 . 8X10"11 m2/s in all three simulations. 123 Chapter Four: RETARDATION OF SORBING SOLUTES IN FRACTURED MEDIA 4.1 Abstract Sorption of a reactive solute during transport in a fractured geologic medium is analyzed under the assumption that within each fracture, retardation varies in proportion to the product of a surface distribution coefficient and the specific surface area of the fracture (i.e. its surface areato-volume ratio). This approach is analogous to the KD model commonly adopted for granular porous media. Numerical migration experiments in discrete fracture networks show that at the plume scale, retardation is both nonuniform and anisotropic. Different segments of the plume, or equivalently different breakthrough fractions at a downstream boundary, are retarded to a different degree. The degree to which various breakthrough fractions are retarded varies as a function of the orientation of the mean hydraulic gradient relative to the orientation of the fracture sets. This variation can be described in the form of a retardation ellipse. The mean arrival time is the only breakthrough fraction that is consistently described by a uniform, isotropic retardation 124 factor. The degree of nonuniformity and anisotropy in the retardation factor is controlled largely by the difference in the mean apertures of each fracture set, the standard deviation in fracture aperture, and the orientation of the fracture sets relative to the mean hydraulic gradient. A larger variability in fracture aperture, or in the orientation of fractures within each set, promotes nonuniform retardation while reducing the degree of anisotropy in the retardation factor. Our results suggest that a transport model based on the conventional advection-dispersion equation, using a uniform retardation factor, may be conceptually incorrect, even for a dense fracture network. While the effects of nonuniform retardation modify the dispersive flux in a way that could be described by "effective" dispersion coefficients, it is our opinion that a model formulation is needed that clearly separates chemical effects from those due to hydrodynamic dispersion. 4.2 Introduction Several countries including Canada, France and Sweden are currently assessing the concept of permanent storage of highlevel nuclear fuel wastes in crystalline rock. The primary pathway for solute transport is likely to be through a network of hydraulically-connected fractures. Any interaction between 125 the dissolved solutes and the fracture walls will result in the retardation of that solute and hence, delay its arrival at the biosphere. Based on laboratory migration experiments of radionuclides in a single fracture, Neretnieks and co-workers proposed two types of retardation mechanisms (Neretnieks et al., 1982); (i) surface sorption at the fracture walls, and (ii) diffusion into and bulk sorption in the rock matrix. Both types of sorption reactions are assumed to be fast, reversible processes which can be modeled using distribution coefficients obtained by equilibrating solutions containing the radionuclides of interest with the rock surface or rock matrix, respectively (e.g. Neretnieks et al., 1982; Vandergraaf et al., 1988). Which of the two sorption mechanisms dominates depends among other factors on the porosity of the matrix, flow rates in the fractures, diffusion coefficients, and the time scale of interest. There is a large potential for bulk sorption due to the large volumes of rock matrix available relative to fracture surfaces (Neretnieks, 1983). However, rates of diffusion into the rock matrix are generally much lower than flow rates within an open fracture. Generally speaking, surface sorption may prevail at small time scales whereas bulk sorption may be dominant at large (geological) time scales. In this chapter we focus on retardation of a solute due to surface sorption on the fracture walls; matrix diffusion and 126 bulk sorption in the rock matrix are ignored for simplicity. The first model of surface retardation was formulated by Burkholder (1976) for the case of a single fault and adapted to a fracture plane by Freeze and Cherry (197 9). According to this model the velocity of a sorbing solute is reduced by a constant retardation factor given by R = l+-^Ka (1) Thus retardation varies in proportion to a surface distribution coefficient, Ka, and the surface area-to-volume ratio of the fracture (SA/V). The determination of the surface area-to-volume ratio of a fracture is difficult in a field setting. A pragmatic approach is to assume two planar, reactive fracture surfaces in an open fracture of aperture b so that the retardation equation becomes (Freeze and Cherry, 1979) * = 1+ f *« (2) Equation 2 indicates that retardation decreases with an increase in fracture aperture. Note that this behavior follows from the surface area-to-volume ratio of a fracture, indicating how much fracture surface area is available for sorption relative to the fluid volume containing the solute. 127 It is not related to the fact that fluid and solutes typically travel faster in wider fractures; the sorption reaction is assumed to be sufficiently fast to be independent of the residence time in a fracture. This sorption model has been used to describe breakthrough curves of sorbing solutes in machined fractures (Vandergraaf et al., 1988), natural fractures in laboratory experiments (e.g. Neretnieks et al., 1982; Vandergraaf and Drew, 19 91), and field experiments (e.g. Abelin et al., 1985). Numerical migration experiments by Moreno et al. (1988) suggested that a similar surface retardation factor could be written for a fracture with variable aperture. They found that for the case of a distributed aperture the mean aperture simply replaces the constant aperture, b, in the formulation of R. To our knowledge, no laboratory or field experiments have yet been performed to evaluate the suitability of an analogous retardation factor at the scale of a fracture network. However, by analogy to the case of a single fracture with variable aperture, one may postulate that the retardation factor in a fracture network is proportional to the average of the surface area-to-volume ratios of all fractures in the network. Henceforth, we refer to this medium property simply as the specific surface area of the fractured medium. Recently, Dverstorp et al. (1992) used numerical migration experiments to study the influence of variable fracture aperture on retardation in sparse fracture networks. The 128 breakthrough curves at a downstream boundary were analyzed with a one-dimensional advection-dispersion equation using a surface retardation factor based on the specific surface area. For a large standard deviation in the aperture distribution, the retardation of the median arrival time was found to be consistently lower than predicted by the retardation model based on a retardation factor using the (known) specific surface area of the network. Based on the observed R the authors calculated effective specific surface areas for the center of mass (median arrival). Dverstorp et al. (1992) concluded that flow channeling significantly reduced the effective values of specific surface area experienced by the solute plume compared to the average value of the fractured medium. However, is the use of a constant retardation factor at the network scale, either based on the average medium value or an effective value of specific surface area, conceptually correct? In other words, can retardation of a sorbing solute plume in a fracture network be considered uniform? First insights can be gained by drawing comparisons to retardation in heterogeneous porous media. Stochastic analyses of reactive solute movement in heterogeneous porous media have shown that spatial variability in the retardation factor may delay the breakthrough, that is, produce a greater effective field-scale retardation factor than expected from any average of the local retardation factors (e.g., Cvetkovic and Shapiro, 1990). 129 Furthermore, negative cross correlation between the local fluid velocity and the local retardation factor may produce enhanced dispersion of the retarded solute plume (e.g., Garabedian et al., 1988/ Valocchi, 1989). Recently, Bellin et al. (1993) developed an analytical solution for the spatial variance of a solute plume in a physically and chemically heterogeneous porous medium. Their results suggest that chemical heterogeneity enhances longitudinal spreading of the retarded solute plume but has little effect on transverse spreading. Numerical migration experiments by Bosma et al. (1993) showed that longitudinal dispersion of a solute plume may be reduced for a positive correlation between sorption coefficient and hydraulic conductivity, relative to the case of negative or no correlation. Fractured media are often extremely heterogeneous with respect to fracture aperture which influences the local retardation factor (see equation (2)). Direct observations on cores (e.g. Snow, 1970) as well as hydraulic measurements (e.g. Raven, 1986) have shown that fracture apertures can span several orders of magnitude in a fractured rock. Since fracture aperture influences retardation at the scale of a single fracture, local retardation factors can also be expected to vary widely in a fractured rock. In addition, there is a conceptual argument for a correlation between the retardation factor and the fluid velocity in a single fracture. The retardation model (2) states that retardation in 130 a fracture is inversely proportional to its aperture. However, fracture aperture also strongly influences the velocity of the fluid. This non-linear coupling of fluid velocity and local retardation may prohibit the use of a constant retardation factor at the network scale. This heterogeneity in local retardation is a direct consequence of the geometric properties of the fracture system; additional chemical heterogeneity will also occur in fractured media but has not been considered here. In this chapter we evaluate the continuum concept of a retardation equation with a constant retardation factor for a fractured medium. First we review the derivation of the continuum transport equation with retardation due to surface sorption and derive a retardation factor based on the specific surface area of a fracture network. Next, we use numerical migration experiments in discrete fracture networks to show that retardation due to surface sorption is a nonuniform and anisotropic process. Different segments of a plume, or equivalently, different breakthrough fractions at a downstream boundary, are retarded to a different degree. Hence surface sorption may either reduce or increase the dispersion of the breakthrough curve depending on the network geometry and its orientation in the flow field. We will show that only the mean arrival time is predicted accurately when using a constant retardation factor based on the specific surface area. 131 4.3 Retardation Models The concept of retardation is based on experimental data which show that the partitioning of many solutes between the solid and the aqueous phase is a fast, reversible process which can be described by a linear isotherm. With the assumption of linear sorption the transport equation for a sorbing solute in a homogeneous porous medium (e.g. Freeze and Cherry, 1979) is given V-(PVC) _ v-(cv) = dc Rf Rf { a) dt * where V is the average linear velocity of the groundwater and D is the dispersion tensor. Rf is known as the retardation factor Rf = 1 +^KD (3b) where 0 and pb are the porosity and the bulk density of the porous medium, respectively, and KD is the distribution coefficient (all assumed to be constant). It is seen from equation (3a) that the effect of sorption is to decrease both the advective flux and the dispersive flux by a constant retardation factor, Rf. The retardation concept may also be applied to a fractured 132 medium where the host rock is assumed impermeable and all sorption reactions are limited to the fracture walls (surface sorption) . Here, the amount of solute, S*, sorbed to the solid in an isotherm experiment is better expressed in terms of the surface area of the sorbate. Introducing the surface distribution coefficient, Ka, we have (Neretnieks et al., 1982) S* = KaC (4) where the mass of solute, S*, sorbed on the solid part of the fractured medium per unit surface area of sorbate is proportional to the solute concentration, C, in solution. Ka can be determined in static sorption experiments by exposing rock coupons containing fracture surface material to a solution containing a sorbing solute of interest (e.g. Vandergraaf et al., 1988). The transport of a reactive solute in a fractured medium of constant porosity can be described in a general form as ,5, -ev.F = * §+ § where F is the total (advective and dispersive) flux of the solute C in a unit volume of rock, 0 is the porosity and q is the mass of solute that is transferred to the sorbate per unit 133 volume of rock. The derivation of the surface retardation factor for a fractured medium is analogous to that for a porous medium with one important exception. We need to consider the bulk surface area, Ab, of the fractured medium in order to relate q in the transport model (5) to S* given by the isotherm experiment (4). In the context of surface sorption in a fractured medium, q is linked to S* via the relation q=AbS* (6) where Ab is defined as the total surface area of all hydraulically-connected (open) fractures per unit volume of rock. Differentiating (4) and (6) with respect to time and assuming 0 and Ab are constant in space and time, the transport equation for a sorbing solute in a fractured medium simplifies to -V.F=(1^,)| (7) The term in parentheses is a constant which we define as the surface retardation factor, <R>, for the fractured medium. The symbol "< >" identifies parameters which are used at the network scale and which need to be distinguished from equivalent parameters at the scale of a single fracture. Assuming all fractures can be approximated by parallel plates 134 of variable length, ln, and aperture, bn, and further assuming a constant depth, D, for all fractures that is much greater than b n (i.e. a two-dimensional problem), the bulk surface area and porosity can be written for a unit volume, V, of fractured rock with N fractures as n 4f Ahb = — with dimensions V ri,2i -»—-± [ L (8) 3] and N T,1rPIP (9) 8 = -x with dimensions V -1= i [L3] We define the ratio of Ab over 0 as the specific surface area, <a>, of the fractured medium N A ^ Q 2 £2« = <a> = —5=1 1 (10) E ^ N A=I We emphasize that <a> is defined as the ratio of the surface area of (open) fractures to the fracture volume, and not to the total volume, V, of the fractured rock. Substituting (10) into (7), the surface retardation factor, <R>, for a fractured 135 medium becomes <R> = l + <a>Ka (11) Note that this continuum retardation equation (11) for a fractured medium is consistent with the retardation equation (2) given earlier for a single fracture. If the fractured medium consists of a single fracture, <a> simply reduces to the surface area-to-volume ratio of this fracture (i.e. 2/b). The specific surface area, <a>, is calculated from all fractures present in a representative elementary volume (REV) of the fractured medium. Dverstorp et al. (1992) have termed <a> the "geometric estimate" of the specific surface area. Similarly, <R> in (11) may also be referred to as the geometric retardation factor. The use of <a> in the continuum retardation equation is justified only if the solute plume encounters a representative (unbiased) sample of fractures along its flow path. Any significant channeling (i.e. a bias towards transport in wider fractures) that involves the entire solute plume should yield a reduction in the observed ("effective") value of specific surface area and consequently the retardation factor. There is considerable practical advantage in replacing the explicit formulation of <a> (equation (10)) by an approximation involving only the mean values of fracture 136 length and aperture as well as fracture density. Consider a network of two fracture sets with N1 and N j fractures, respectively. Assuming that the fracture aperture and fracture length of each fracture set are independent random variables given by a probability distribution with respective mean values 1J-, b / for fracture set I, and lm3, bm3 for fracture set J", respectively, the specific surface area can be approximated as <a> = 2N±li ~ N+ + 2Njli ~ (12) Vb^+NU'b" It should be kept in mind that the independence of the fracture length and aperture is a crucial assumption in identifying this approximation for <a>. This assumption may or may not be satisfied in fractured geologic material depending on the genesis of the fracture network (e.g. Chernyshev and Dearman, 1991). Equation (12) is easily rewritten in terms of fracture density. The fracture density, y, of a fracture set can be expressed by the ratio of the total number of fractures in this set per volume of rock (e.g. yi=Ni/V) . Substituting the product of the respective fracture density and V for total numbers of fractures of each set in (12), V cancels and we have 137 2 Y i l m i + 2 Y J l,n J ; : •—; r <a> = (13) This formulation of specific surface area <a> is readilyextended to three dimensions by assuming a probability distribution for fracture depth (instead of a constant D ) . This extension would introduce the mean fracture depth of each fracture set into equation (13) . In the following we consider a few special cases for which the specific surface area <a> of a network with two fracture sets can be further simplified. These special cases will be explored in our migration experiments to be discussed later. If the two fracture sets have a common mean aperture bm, that is b,,,1 is equal to bmi, we have <a> = — V Y m . Y mJ . = ~ (14) For this special case, <a> can be estimated directly from the common mean aperture. The specific surface area is independent of the mean fracture length and the total number of fractures of the two fracture sets, respectively. In contrast, if we assume that the two fracture sets have a common mean fracture length lra, that is lml is equal to lmJ, we have 138 2I <a> = ( Y i+ m5L1I ; Y:') , = —=-A!— 2 ( v i + YY:') , (15) Here <a> is given by twice the inverse average of the two mean aperture values weighted by the relative fracture densities of the two fracture sets. For a fracture network composed of two fracture sets with equal fracture density and a common mean fracture length lm, the specific surface area is given by <a> (16) that is twice the inverse average of the two aperture means. In summary, the transport equation for a sorbing solute in a fractured medium is given as V-F = (l+<a>JCa) J£ (17) where the term in parentheses is defined as the (geometric) surface retardation factor, <R>. For a fracture network composed of two fracture sets the specific surface area is given by the general formulation (13). The above transport 139 equation applies to a fractured medium that satisfies the following assumptions; (i) fracture planes are open and can be approximated by parallel plates, (ii) fracture length and fracture aperture are statistically independent, (iii) no significant channeling (biased sampling of wide fractures) of the entire solute plume occurs at the scale of the REV, and (iv) a uniform and isotropic retardation factor exists at the network scale. The last assumption is perhaps the most interesting. It implies that at the continuum scale the entire flux of the solute (advection and dispersion) has to be reduced equally by sorption. In addition this retardation of the solute plume has to be independent of the direction of the hydraulic gradient. Only then is the use of a constant retardation factor <R> justified. Most of the discussion in this chapter will focus on assumption (iv). A transport experiment in the field or laboratory is needed to test the validity of these assumptions. For this purpose the retardation factor would be determined experimentally using the retardation equation for a fractured medium <R> = ^ ^ = l + <a>JC_ a <vr> v (18) ' where <v> and <vr> are the average linear velocities of the groundwater (or a conservative solute) and a reactive solute, respectively, in the fractured medium. Commonly, retardation 140 factors are determined from breakthrough curves (BTC) in tracer experiments. Brusseau and Rao (1989) pointed out that the mean arrival time is the correct reference point to estimate a retardation factor. This method yields good estimates of <R> even in the case of non-ideal (nonequilibrium) sorption behavior because the first time moment is independent of enhanced spreading due to kinetic limitations in the sorption behavior (Brusseau and Rao, 1989 and references therein). Although (18) provides a way to estimate <R>, it does not test whether the same retardation applies to all segments of the solute plume travelling through the medium. Assuming that molecular diffusion does not contribute significantly to dispersion, a more general retardation equation can be written <v>rof <R> ref ref = — <v > * v where the index ref r = l + <3>X a (19) a ref indicates any reference point of the solute plume (or equivalently any reference point on a BTC at a downstream boundary). The general retardation equation (19) above is satisfied if the entire breakthrough curve is retarded equally, that is shifted by a constant <R>. This relation implies that surface sorption reduces the advective and dispersive flux by a constant <R> as assumed for the 141 reactive transport equation (17). The numerical migration experiments described below will demonstrate that this assumption is not generally justified. 4.4 Numerical Migration Experiments: Design and Choice of Parameters The fractured medium in our migration experiments is represented by two-dimensional fracture networks generated stochastically from a set of fracture system statistics (Clemo 1994). The numerical model generates a network of open fractures by randomly locating fracture midpoints using a uniform distribution. Fracture orientations are normally distributed and fracture length is specified by a negative exponential probability distribution. Fracture aperture is log-normally distributed and does not vary along the length of a fracture. These stochastic input parameters are all independent. The networks studied here consist of two sets of fractures. Each set is specified by scan-line density (ys) , mean length (lm) , mean aperture (bm) and the standard deviation (a) of aperture in log10 space. Scan-line density is the expected number of intersections of the fractures in a set with a line of unit length that is normal to the mean orientation of the set. Flow and transport is simulated using a modified version of 142 the finite difference model DISCRETE (Clemo, 1994). The fracture segments between fracture intersections have a constant flow. The flow is related to the hydraulic head drop between intersections by 12 JJ Al where Q is the volumetric flow in a fracture segment of unit depth, Ah is the head drop between the intersections, and Al is the distance between the intersections. Conservation of mass at the intersections results in a set of linear equations that describe the flow within the domain. Transport is modeled using a particle tracking routine. For any single network realization, two runs are performed; (1) transport of a non-sorbing ("conservative") solute and (2) transport of a reactive solute sorbing to the fracture walls (equilibrium surface sorption). Non-reactive particles travel from intersection to intersection according to the residence time of the fluid, x, in a fracture segment given by r = b-Al/Q (21) Streamtube routing is implemented at the fracture intersections. Surface sorption is modeled by increasing the 143 residence times of all fracture segments by their respective surface retardation factor, R, as defined by the retardation equation for a single fracture (equation (2)). This transport routine assumes plug flow, that is, dispersion at the single fracture scale is negligible compared to dispersion at the network scale (e.g. Hull et al., 1987). This assumption is a reasonable approximation for the vast majority of fracture segments under the flow conditions studied in our migration experiments. The flow domain and the boundary conditions are chosen to mimic the conditions of a tracer experiment at the near-field scale. The flow domain is 10m by 10m with constant head boundaries on the left and right and no-flow boundaries at the top and bottom (Fig. 4.1). The magnitude of the head gradient between the constant head boundaries was set to 0.01. However, the magnitude of the global hydraulic head gradient does not influence the retardation behavior studied here. Several migration experiments in dense networks (see reference geometry below) were rerun on a larger flow domain (20m by 20m) with identical boundary conditions as those shown in Figure 4.1. No significant differences in the retardation behavior were observed, suggesting that a domain size of 10m by 10m is sufficiently large to study a representative retardation behavior in the dense networks considered here. The solute is introduced as an instantaneous pulse along a line source of one meter length in the center of the left (up144 gradient) boundary (Fig. 4.1). The line source is centered in the middle of the inflow boundary so as to minimize the influence of the no-flow boundaries on the spreading of the solute plume. A source length of one meter was sufficiently large to provide good access (-4-6 injection nodes) to the dense fracture networks studied here. Transport calculations are not performed in individual network realizations if (i) no fracture intersects the source area or (ii) the combined flow rate for all intersecting fractures is less than 10 percent of the flow rate across the domain. This filter ensures that the particles are not trapped in close proximity to the source due to poor connection of the source line to the network. Sensitivity analyses indicated that, although this filter artificially reduces the variance in the transport results (in particular for the sparse network in scenario 2) , it does not influence the general pattern of behaviour we observed. Transport of a single particle is completed when it reaches the down-stream (right) boundary. Its total transport time is given as the sum of the fracture residence times determined during transport through the network. The distribution of total transport times for a large number of particles constitutes the breakthrough curve (BTC) for a given transport simulation. The particle tracking technique requires the use of a sufficiently large number of particles so that the transport characteristics do not change with a further increase in the 145 number of particles. In the networks studied here 5,000 particles are sufficient to characterize the reactive transport behavior. No significant deviations in transport results were obtained in comparative runs with up to 50,000 particles. Four retardation factors, <R>5, <R>50/ <R> m^ anci <R>95 are determined for each network realization by comparing the arrival times of four reference mass breakthrough fractions, T5, T50, Tm, and T95 in a conservative and reactive transport run. These reference points (5%, 50%, mean, and 95% mass breakthrough fraction) have been chosen to compare retardation of the leading edge, the center, and the tail of the solute plume. Note that the solute breakthrough curves in most fractured media are skewed to the right, so that the mean arrival time is later than the arrival time of the center of mass (i.e. Tm > T50) . In accord with the general retardation equation (19) a continuum transport equation with a constant retardation factor can be formulated if these four <R>ref are equal (assuming no diffusion into the rock matrix). In this case of uniform retardation, the four <R>ref may either approximate the (known) geometric retardation factor <R> (equation 11) or converge to an (unknown) effective retardation factor. However, if the four <R>ref differ significantly, the retardation process is considered nonuniform and a constant retardation factor is not applicable. An alternative interpretation of the experiments would be to fit an 146 advection-dispersion model to the BTC of the sorbing solute (e.g. Dverstorp et al., 1992). In this procedure, nonuniform retardation could potentially be accounted for by fitting an effective dispersion coefficient to a retardation run. As discussed in a later section, from a conceptual viewpoint, we prefer to separate hydrodynamic dispersion from that due to chemical interaction. Transport properties (non-reactive and reactive) vary considerably among individual realizations of a given fracture network description. This variability results from differences in the arrangement and connectivity of fractures that are possible within the stochastic network description. Trial runs with up to 500 Monte Carlo runs showed that stable average values for Tref and <R>ref were obtained for -100 Monte Carlo realizations of a given fracture network description. All results presented here are based on 250 Monte Carlo realizations. The extent of retardation due to surface sorption depends on the type of network studied. In this chapter, we study primarily dense fracture networks in which a continuum approximation to reactive transport would appear to be most promising. Figure 4.2 shows an example of a fracture network with the reference geometry used in the majority of the migration experiments. This reference geometry consists of two orthogonal fracture sets each with a mean length of one meter and a scan-line density of six fractures per meter. This 147 c o m b i n a t i o n o f fracture length a n d s c a n - l i n e density r e p r e s e n t s w e l l - c o n n e c t e d a n d dense n e t w o r k s (10-20 intersections/m2) . F r a c t u r e l e n g t h , fracture density a n d fracture a p e r t u r e all i n f l u e n c e s p e c i f i c surface area presumably retardation. (see e q u a t i o n (13)) a n d H o w e v e r , in this chapter w e consider p r i m a r i l y t h e influence of fracture a p e r t u r e on r e t a r d a t i o n . T h e r e f o r e , all n e t w o r k s studied here a r e c o m p o s e d of t w o f r a c t u r e sets w i t h equal m e a n fracture length, l m , of o n e m e t e r a n d equal fracture d e n s i t y . In this way, specific surface area a n d t h e g e o m e t r i c r e t a r d a t i o n factor c a n b e e s t i m a t e d d i r e c t l y from t h e m e a n a p e r t u r e of t h e t w o sets (equation (16) ) . The range of f r a c t u r e a p e r t u r e s , b , in e a c h fracture set is r e p r e s e n t e d b y a lognormal density d i s t r i b u t i o n w h i c h is given (for b > 0 ) as ,, > n(b) = 1 bj2n(lnl0a)2 -(logij-logbj2 — exp-^—„ 2 ,«,„. (22) 2a and specified by two parameters, the mean (log b0) and standard deviation (a) of aperture in log-10 space. The nontransformed value b 0 is the most probable aperture value (geometric mean) and is related to the (arithmetic) mean aperture, b m , by 148 , , (ln(10)a)2 / O ox As can be seen from (23) , the geometric mean, b 0/ is always smaller than the arithmetic mean, b m . Note that the specific surface area (equations (13)-(16)) is estimated based on the arithmetic mean aperture, bm. We study the sensitivity of retardation to the spread (or variability) in the aperture distribution by varying c However, the resulting aperture distributions differ considerably in shape depending on whether the geometric mean b 0 or the arithmetic mean bm is held constant (see e.g., Tsang et al., 1988) . The use of a constant bm generally gives more weight to the lower aperture values relative to that using a constant b 0 . Note that for o=0 the aperture distribution collapses to a single aperture value (b=b0=bm) . In our sensitivity studies the arithmetic mean aperture, b m/ is held constant so that the specific surface area and the geometric retardation factor of the network remain constant. Results of migration experiments for six types of fracture networks (scenarios) are presented in this chapter. The different network parameters used for these scenarios are summarized in Table 4.1. Scenarios 1 and 3-4 use the dense and well-connected reference geometry of two orthogonal sets shown in Figure 4.1. The only parameters to be varied relate to the aperture distribution for the two sets. In the first 149 scenario we study a fracture network where the two fracture sets are described by a common aperture distribution with a constant mean bm equal to 80 jim (Table 4.1) . In scenarios 3-4, we study a network where each fracture set has a distinct aperture distribution with bm equal to 2 0 pm and 80 Jim, respectively. Scenarios 3 and 4 differ in the orientation of the network with respect to the overall direction of flow (Table 4.1). In the remaining three scenarios we study "connectivity effects" by changing the network geometry itself. In scenario 2 we reduce the fracture density of the fracture network with common aperture distribution by half relative to the reference geometry. In the last two scenarios we study fracture networks in which the two fracture sets (with differing mean apertures) are not strictly orthogonal (scenario 5) or show no preferred mean orientation (scenario 6) . For each scenario, arrival times and <R>ref are calculated repeatedly while systematically increasing o from 0 to 0.5 in both sets. The sorption strength of the surface reaction characterized by Ka is held constant in all simulations. We have chosen Ka somewhat arbitrarily (Ka= 4-10"3m) so that the geometric retardation factor <R> is equal to 101 in the first scenario (with <a>=2/bm=2 .5-10+4m) . The chosen value of Ka is within the range of Ka values measured in static sorption experiments for weakly sorbing radionuclides (e.g. Vandergraaf and Drew, 1991) . 150 4.5 Results and Discussion 4.5.1 Two Orthogonal Fracture Sets with a Common Aperture Distribution In the first scenario we study the retardation behavior in a dense fracture network where the two fracture sets are described by a common aperture distribution. Figure 4.3 summarizes the transport behavior of a conservative and sorbing solute in scenario 1. The upper panel shows the average arrival times, Tref, of the conservative solute at the down-gradient (right) boundary as a function of a. The respective retardation factors, <R>ref/ panel. are shown in the lower The arrival times of the sorbing solute are easily back-calculated by multiplying <R>ref by Tref. All values are average values of 250 Monte Carlo realizations with a standard deviation indicated as one-sided error bars. Two general observations can be made in this experiment with respect to retardation at the plume scale. First, the four <R>ref are not equal except for the special case where all fractures have the same aperture (G=0). This behavior suggests that the concept of a single retardation factor rigorously applies only to a network composed of fractures with equal aperture. Provided there is any variability in fracture aperture (<J>0) retardation at the plume scale cannot be 151 considered uniform. The difference in retardation factors <R>ref can be significant depending on the spread in the aperture distribution. For example, consider the fracture network in which fracture aperture varies by approximately one order of magnitude (<3=0.3). In this case the 5 percentile breakthrough fraction is retarded by a factor of -70 whereas the 95 percentile breakthrough fraction is retarded by a factor of -120 (Fig. 4.3). Second, the geometric estimate of the retardation factor (here <R>=101) which is based on the specific surface area, <a>, accurately predicts the retardation of the mean arrival time for the entire range of G values considered (Fig. 4.3). At the same time, the geometric estimate of <R> consistently overestimates retardation for breakthrough fractions ahead of the mean arrival time and underestimates retardation for breakthrough fractions arriving later than the mean arrival time (Fig. 4.3). This systematic deviation of the observed <R>ref from the geometric <R> is a result of biased "sampling" of fractures by the solute particles. For example, a fast particle representing the 5 percentile breakthrough fraction preferentially travels in fractures of wide aperture since these fractures tend to have higher flow rates (see equation (20)) and shorter residence times. It follows that the effective specific surface area experienced by a fast particle is lower than the geometric estimate of <a>; hence the resulting retardation factor <R>5 is also lower than the 152 geometric <R>. Similar arguments can be made to explain the greater retardation of slow particles (e.g. <R>95) relative to the geometric estimate. The mean arrival time Tm of the conservative solute is greater than the arrival time of the 50% breakthrough fraction T50 indicating a BTC which is skewed to the right (Fig. 4.3) . Consistent with the above line of reasoning, the faster center of mass shows a lower <R>50 than the mean arrival time (<R>m~<R>=101) for all o values greater than -0.1 (Fig. 4.3). Such differences cannot be confirmed in a statistical framework due to the high standard deviations in the mean values of Tref as well as <R>ref in particular for large values of a. However, the trend of a lower <R>50 relative to <R>m has been observed in most realizations of this scenario. The observed differences in the values of <R>50 and <R>m raise a question concerning which reference point is more appropriate for estimating the "average" retardation of a solute plume. Brusseau and Rao (1989) favor the use of the mean arrival time to estimate a retardation factor. They argue that Tm is a truer measure of retardation than is T50 because the retardation of Tm is independent of any non-equilibrium effects (kinetic limitations in the sorption reaction). As shown here, retardation of Tm is independent of any spatial heterogeneity of local retardation (nonuniform sorption) in a network of fractures (Fig. 4.3). In principle this allows a prediction of retardation of the mean arrival time from the 153 specific surface area of the fractured medium (or vice versa). From a practical viewpoint, however, this reference time is influenced significantly by the tail of the BTC which can be very long and difficult to measure accurately. This dependence on the tail of the BTC makes it much more difficult to determine Tm in an experiment than the arrival time of the center of mass (T50) . More importantly, the strong weight that the tail of the plume gives to Tm (and subsequently <R>m) results in an estimate of plume retardation which is often not representative of the retardation of the faster "bulk" of the plume. To illustrate this point we have plotted snapshots of the spatial distribution of the particles in a network realization of scenario 1 for the conservative run and the retardation run (Fig. 4.4). The spatial distribution of the conservative solute plume is shown for the time T0 equal to 2 days after pulse injection. At this time all particles are still present in the domain. This spatial distribution should be approximately reproduced by the retardation run at a time t0 multiplied by the "average" <R>. A visual comparison of the retarded plumes for times equal to 154 days (T1=To*<R>50) and 202 days (T2=T0*<R>m) with the conservative plume shows that <R>50 is a better estimator for the retardation of the "bulk" of the plume (Fig. 4.4). In particular, the retardation of the front of the plume is much better captured using <R>50. In contrast, the differences in using <R>50 or <R>m for predicting 154 the retardation of the tail of the plume are minor. In densely fractured media, spatial moments provide a meaningful description of the spatial distribution of the solute plume, that is the first and second moment represent the center of mass and the degree of spatial dispersion of the solute plume, respectively (Smith et al., 1990) . A comparison of the first and second spatial moments of the three solute plumes (Ml and M2, respectively; Fig. 4.4) confirms the visual impression that a retardation factor estimated from T50 describes the spatial distribution of a reactive solute plume (both center of mass and spreading) much better than does an estimate of <R> based on Tm. Unfortunately, the geometric estimate of <R> is not a good model for <R>50, in particular for high a (Fig. 4.3). 4.5.2 A Continuum Model for Nonuniform Retardation? The results of scenario 1 clearly demonstrate that even for a dense and well-connected fracture network retardation is more complicated than the conventional retardation model applied in a granular porous medium. The major difference between- the Ka approach in a fractured medium and the KD approach in a porous medium lies in the way local retardation is linked to local particle velocities. This point is best understood by considering the influence of sorption on the 155 particle velocity distribution. In the Fickian dispersion model for a porous medium, the solutes are assumed to sample randomly from a velocity distribution of which the first and second moment define the average linear velocity and the dispersion of the conservative solute, respectively (Bear, 1972). In the traditional KD approach, it is assumed that the entire velocity distribution is shifted by a constant retardation factor for a sorbing solute. In this way, the average linear velocity and the dispersion coefficient are both reduced by this retardation factor (equation 3a). Because it is not possible to measure the velocity distribution this assumption can only be tested at the continuum scale using an equivalent form of the general retardation equation (19). In the Ka approach developed here, we assume a priori a relationship between the velocity of a conservative solute and a sorbing solute at the scale of a single fracture in the form of equation (2). In principle, the reactive particle velocity distribution is shifted by a constant factor only if aperture is constant for all fractures of the network, that is O is equal to 0. The use of discrete fracture networks allows us to calculate the distribution of local particle velocities in the network. Figure 4.5 shows such particle velocity distributions for a conservative and a retardation run observed in a single network realization of scenario 1. As expected, for a equal to zero the entire velocity distribution for reactive particles is shifted towards a lower average particle velocity but the 156 spread in particle velocities is constant (Fig. 4.5a). In contrast, for a fracture network with variable aperture (e.g. 0=0.5) the velocity distribution is not simply shifted by a constant factor; instead the mean as well as the standard deviation of the reactive particle velocity distribution is changed relative to the conservative particle (or fluid) velocity distribution (Fig. 4.5b). This perturbation in the velocity distribution arises because fracture aperture influences not only the local retardation factor but also the fluid velocity itself (see equation (20)). In scenario 1, this coupling of fluid velocity and retardation via fracture aperture gives rise to a substantial increase in the spread of the velocity distribution when introducing retardation (Fig. 4.5b). This increased spread in the particle velocities explains the enhanced dispersion of the retarded solute plume relative to the plume of a non-sorbing solute (Fig. 4.4). Dispersion of the retarded solute plume remains unchanged only in a network with no spatial variability in local retardation (uniform retardation for o=0.0, Fig. 4.3) in which the spread in the particle velocity distribution is unchanged (Fig. 4.5a). This enhanced dispersion in the retarded solute plume due to a coupling of local retardation and local fluid velocities is not limited to a fractured medium. Stochastic analyses of reactive solute movement in a heterogeneous porous medium have shown that a negative cross correlation between the pore water 157 velocity and the local retardation factor may produce enhanced dispersion of the retarded solute plume (e.g., Garabedian et al., 1988/ Valocchi, 1989/ Bellin et al., 1993/ Bosma et al., 1993). It should be kept in mind, however, that our results are not directly comparable to the above studies which all assume a correlation between the two random fields of sorption coefficients (or retardation factors) and hydraulic conductivity. In the Ka approach for a fractured medium, the fluid velocity and retardation factor in a single fracture are related via the fracture aperture. Yet, the overall nature of this relationship is more complicated than a simple negative cross correlation. The fluid velocities are also influenced by the local head gradients which in turn depend on the network geometry as well as the boundary conditions. Hence the solution to the flow problem has a direct bearing on the nature of the relationship between fluid velocity and retardation. Bellin et al. (1993) further suggested that chemical heterogeneity has a much greater effect on longitudinal spreading than on transverse spreading of the retarded solute plume in a heterogeneous porous medium. A preliminary analysis of the spatial moments in scenarios 1-2 suggests that this trend may also exist in a fractured medium (e.g. compare spatial moments in Figure 4.4). A more detailed analysis is required to confirm this hypothesis. In this context it needs to be emphasized that the spatial heterogeneity in retardation 158 in our fracture networks is a result of the physical heterogeneity (i.e. variable aperture) and not chemical heterogeneity (that is variable KD or Ka) as assumed by Bellin et al. (1993) . The presence of a variable sorption strength (Ka) has not been considered here but would certainly further complicate the process of nonuniform retardation. How can we describe this nonuniform retardation at the continuum scale? The non-linear coupling of retardation and fluid velocity prohibits in all but a few special cases the use of a transport equation with a constant retardation factor as postulated in equation (17). However, to our knowledge a rigorous mathematical model that describes the variable retardation observed (e.g. Fig. 4.3 and 4.4) in a transport equation is currently not available. A pragmatic alternative would be to consider only the retardation of the "advective flux", that is retardation of the first time moment (or perhaps retardation of the first spatial moment) of the solute plume. The effects of nonuniform retardation for non-central parts of the solute plume would then be absorbed in the operator describing the "dispersive flux" using a lumped parameter. This approach is essentially pursued when fitting the advection-dispersion model (equation (3)) to breakthrough curves of a sorbing solute without obtaining the dispersion coefficient a priori from a conservative transport experiment (see Dverstorp et al., 1992 for such an application to a fracture network). For example, in scenario 1 with o equal to 159 0.5, nonuniform retardation results in enhanced dispersion and a fit of the advection-dispersion equation would yield a greater "effective" dispersion coefficient for the sorbing solute than for the conservative solute. Although appealing in its simplicity, we do not recommend this modelling approach for conceptual reasons. The use of a single, lumped parameter which would describe the physical as well as chemical processes (sorption) resulting in dispersion introduces ambiguity into the meaning of the dispersion coefficient. The dispersion coefficient would no longer be solely determined by the hydromechanics of the medium (such as dispersivity) but it would also depend on the sorption strength of any given solute. From a pragmatic point-of-view this would imply that tracer tests would have to be run for every sorbing solute of interest to determine its "effective" dispersion coefficient. In our opinion, it is preferable to find a model formulation for reactive transport in which each parameter is defined as much as possible by a unique physical or chemical process. Our analysis of the influence of fracture network parameters on nonuniform retardation, discussed in the remainder of this chapter, may help to develop a continuum model that accounts explicitly for nonuniform retardation. 160 4.5.3 Influence of Fracture Density on Retardation In scenario 2 we have reduced the scanline density of both fracture sets to 3 fractures/m, that is half the scanline density of scenario 1 (Table 4.1) to demonstrate the influence of fracture density on nonuniform retardation. Figure 4.6 summarizes the transport behavior in this relatively sparse network type. As a result of the reduction in the number of fractures present in the domain, the conservative travel times are generally much greater than in scenario 1 (Fig. 4.6a). Apart from this expected reduction in the permeability of the fracture system, the overall transport behavior is remarkably similar to that described for scenario 1. Again we observe a growing spread in the conservative arrival times Tref for larger a (Fig. 4.6a). The retardation of the mean arrival time is predicted by the geometric estimate of <R> (Fig. 4.6b). Note, that the value of <R> is the same for scenarios 1 and 2 (<R>=101) despite the difference in total numbers of fractures. This behavior follows since <R> is governed by the specific surface area <a> which should be independent of the total (or relative) number of fractures present in the networks of scenarios 1-2 which have a common mean aperture (where b m 1= b m 2 , see (14) ) . As observed earlier, the various <R>ref deviate from the geometric retardation factor for o greater than zero (Fig. 4.6b). Despite these similarities, the transport behavior in 161 scenarios 1 and 2 differ in two important aspects. First, a reduction in fracture density introduces a larger variability in the estimates for <R>ref and particularly for the conservative Tref, as indicated by the larger errorbars for the average values of <R>ref and Tref (Fig. 4.6) . Second, a sparser network reduces the degree of nonuniformity in retardation, especially at larger o (compare Fig. 4.3 and Fig. 4.6). Both observations are related to the extent of flow channeling that occurs in the respective network geometries. An increase in o generally provides more opportunity for flow channeling due to the increased variability in fracture aperture (e.g. Nordqvist et al., 1992/ Dverstorp et al., 1992). However, the continuity and dominance of individual flow channels differ for different fracture densities. Figure 4.7 shows the flow distribution for representative single realizations of scenario 1 (top) and scenario 2 (bottom) for o equal to 0.5. The relative thickness of the bar indicates the amount of flow (volumetric flow rate) passing through a hydraulically-connected fracture segment. The flow distribution in these two realizations is typical in that the dense and well-connected fracture network develops several discontinuous flow channels throughout the domain (Fig. 4.7a). In contrast, the sparse and poorly connected fracture network develops one or two dominant and rather continuous flow channels across the flow domain (Fig. 4.7b). This difference in the type of flow channeling has 162 implications for both conservative and reactive transport. In a sparse network the extent of a dominant flow channel and its exact position with respect to the source area is crucial for the conservative transport times. Good access of the source area to the dominant channel will result in fast travel times with little dispersion of the conservative BTC; poor access will result in long travel times and a conservative BTC with a very long tail. When averaging over many Monte Carlo realizations, this uncertainty in the exact position of the dominant flow channel translates into large standard deviations of the average values of Tref (Fig. 4.6a) . Reactive transport is affected similarly by the dominant flow channels developing in a sparse network. Transport of sorbing solutes in such a channel results in less opportunity for biased sampling of fracture apertures. Slow particles as well as fast particles will use this flow path for much of their travel across the domain. As a result, at the scale of individual fractures, they sample a similar suite of surface area-to-volume ratios and thus the retardation factors of the various reference mass fractions are similar as well. Typically, a good connection of the source area to a dominant flow channel will result in low <R>ref with little nonuniformity whereas a poor connection to a dominant flow channel (or the absence of the continuous flow channel) will give rise to a greater spread of <R>ref. The overall effect on retardation at the plume scale is a reduction in the spread of 163 <R>ref albeit with higher variability in the average values of <R>ref (Fig. 4. 6b) . In contrast, in a dense network with discontinuous flow channels, the solute plume disperses more widely into the network, sampling many more fractures. The discontinuous flow channels are only of local importance in focussing flow and transport. The resulting wider spreading of the solute plume and sampling of many more fractures (with different surface area-to-volume ratios) results in a higher degree of nonuniformity in retardation at the plume scale. In a dense network, the exact position of the (discontinuous) flow channels with respect to the source area is not crucial and estimates of <R>ref and Tref are less variable. We emphasize again that the reduction in fracture density and the resulting tendency for extended flow channels did not cause a general lowering in retardation of the plume. This result appears at first to contradict findings by Dverstorp et al. (1992) obtained in sparse fracture networks using numerical migration experiments. They observed that retardation of the center of mass (T50) was on average a factor of two lower than expected based on the geometric estimate of specific surface area. From this, the authors concluded that flow channeling reduces the specific surface area. However, Dverstorp et al. failed to acknowledge the nonuniform character of retardation. It is conceivable that the observed reduction in retardation is due entirely to nonuniform 164 retardation, that is, a lower retardation for the center of mass relative to the mean arrival time. According to our results significant channeling of flow and transport does not, averaged over many Monte Carlo realizations, reduce the retardation of any segment of the solute plume. Instead, the major effect of strong channeling is to increase the variability in <R>ref for a given network description. The large variability observed in the (conservative and reactive) transport behavior in the sparse network suggests that the domain size of 10m by 10m is below the scale of a representative elementary volume (REV) for such a low fracture density (if this REV exists at all). We have not studied larger domain sizes for scenario 2. Nevertheless, it is reasonable to conclude that a decrease in fracture density generally reduces the degree of nonuniformity in retardation. However, there is a greater uncertainty attached in estimating retardation of any segment of the plume. 4.5.4 Two Orthogonal Fracture Sets with Differing Mean Aperture In scenarios 3-4 we study retardation in an orthogonal fracture network consisting of one fracture set with wide fractures (bm=80]Lim) and a second set with relatively narrow fractures (bm=20jim) . Here the orientation of the network with 165 respect to the mean hydraulic gradient has an effect on retardation. Figure 4.8 shows the retardation behavior for the two scenarios. Note that the only difference between scenarios 3 and 4 is a change in orientation of the entire network by 90 degrees within the flow domain (Table 4.1). Three observations summarize the retardation response in a network where the fracture sets do not have a common aperture distribution. First, the observed retardation factors <R>ref in both scenarios are not equal for the case of two constant aperture values in the respective sets (0=0/ Fig. 4.8). This behavior confirms our earlier contention that retardation is only uniform in a network where all fractures have the same aperture value. Second, retardation of the mean arrival time, <R>m, is again described accurately by the geometric retardation factor for both orientations of the fracture network (Fig. 4.8). The geometric estimate of <R> in this network (here <R>=161) is greater than in scenario 1 (<R>=101) because the specific surface area in scenarios 3-4 is greater due to the smaller mean aperture of the second set (see equation (16)). By definition, the geometric <R> is independent of the orientation of the mean hydraulic gradient to the network and hence is equal in scenarios 3 and 4. The good agreement of this geometric <R> and <R>m demonstrates that for these two orientations of the hydraulic gradient, the retardation of the mean arrival time is controlled by the specific surface area 166 of the network. We will show below that <R>m agrees with the geometric <R> for any orientation of the fracture network relative to the mean hydraulic gradient. Third, the difference in the various <R>ref for a given scenario is not only a function of G but also depends on the orientation of the network relative to the mean hydraulic gradient. In scenario 3, where the fracture set with the larger bm is oriented in the direction of the mean hydraulic gradient, the degree of nonuniformity increases as a function of G (Fig. 4.8a). The geometric <R> systematically overestimates <R>ref for reference breakthrough fractions ahead of Tm and underestimates <R>ref for reference breakthrough fractions arriving later than Tm (Fig. 4.8a). This behavior is a result of preferential travel of fast particles in wide fractures as discussed for the case of scenario 1. Here the nonuniformity in retardation is more pronounced since the range of fracture apertures is greater (note the change in the scale of the <R>ref axis between Figures 4.3b and 4.8) . In scenario 4, when the fracture set with the smaller bm is aligned with the direction of the mean hydraulic gradient, the observed <R>ref fall above and below the geometric <R> when considered over the entire range of G (Fig. 4.8b). For small G, the geometric <R> underestimates the <R>ref for reference breakthrough fractions ahead of Tref and overestimates <R>ref for reference breakthrough fractions behind Tref (Fig. 4.8b). This pattern is gradually reversed as G increases. 167 The value of C equal to ~0.2 represents a cross-over point at which the four <R>ref are equal. But what process causes the early breakthrough fractions to be retarded more than the late breakthrough fractions for small c in scenario 4? This observation is not intuitive and is "forced" by the orientation of the two fracture sets with respect to the mean head gradient. Consider the case of two fracture sets with constant but differing apertures (b1=20inm and b2=80p.m for c=0) . In this configuration, fluid velocities are often greater in the narrow (horizontal) fractures as the influence of the small cross-sectional area (aperture) is more than -compensated for by high local head gradients. Similarly, fluid velocities are often smaller in the wider (vertical) fractures due to the small local head gradients. As a result fast particles will have spent the majority of their (relatively direct) travel path across the domain in horizontal fractures with high retardation. In contrast, slow particles use a more circuituous path involving many more vertical fractures and thus experiencing relatively less retardation. With an increase in the variability of fracture aperture in both sets (larger a) the difference between wide fractures in one set and narrow fractures in the other set becomes less distinct. At some point (here O£0.2) there is a sufficient number of wide fractures in the horizontal set and narrow fractures in the vertical set so that the fastest particles again travel mostly in wide, horizontal fractures 168 experiencing the lowest retardation. Note that the unusual form of nonuniform retardation (earlybreakthrough fractions are retarded more than the late breakthrough fractions) for small (J in scenario 4 implies a reduced dispersion of the retarded solute plume relative to a conservative solute plume. Such a reduction in plume spreading has also been observed in a heterogeneous porous medium for the case of a positivily correlated sorption coefficient KD and hydraulic conductivity (Bosma et al., 1993) . The reason for this reduced dispersion is essentially the same in both fractured and (heterogeneous) porous media, namely a positive correlation of the local retardation factor and the local fluid velocity. However, the underlying cause for this correlation is physical heterogeneity in the case of the fracture network (variable aperture) and chemical heterogeneity (variable KD) in the case of a heterogeneous porous medium. This leads to different retardation behavior in response to a change in the physical heterogeneity of each system. For example, an increase in aperture variability (increase in a) in scenario 4 produces once again enhanced dispersion of the retarded solute plume indicative of a negative correlation between local retardation factors and local fluid velocities. This change in correlation structure between local retardation and fluid velocity as a function of physical heterogeneity is perhaps unique to a fractured medium. 169 The different retardation response in scenarios 3 and 4 demonstrates the directional dependence of retardation for various segments of the plume. Therefore, retardation at the plume scale has to be considered anisotropic in networks in which the mean fracture aperture differs between two (or more) fracture sets. Given the orientation of the network in these two scenarios (0° and 90° angle from the direction of the mean hydraulic gradient), one might expect that these results represent the maximum possible range of <R>ref due to rotation of the network alone. In order to test this hypothesis, two fracture networks (G=0.0 and 0.5) of the same network geometry as in scenarios 3 and 4 were rotated in increments of 15 degrees through the flow field using identical boundary conditions as shown in Figure 4.1. For each orientation, 250 Monte Carlo realizations were used to determine average values of <R>ref. The four <R>ref determined in these simulations are plotted in polar coordinates, where the angle represents the degrees of rotation of the network and the radius is given by the respective <R>ref (Fig. 4.9). Figure 4.9a shows results for o=0.0, while Figure 4.9b is for a=0.5. In this presentation the four <R>ref for scenario 3 (^=0°) plot on the x-axis and the <R>ref for scenario 4 (<x1=90°) plot on the y-axis. These two orientations are referred to as the principal orientations of the network with respect to retardation. Similarly, the <R>ref measured in those two orientations are defined as the 170 principal retardation factors for the respective reference point. The simulation results show that these principal retardation factors indeed represent the maximum range for any given <R>ref (Fig. 4.9) . The <R>ref observed for the intermediate orientations are predicted well by a ellipse retardation using the observed principal retardation factors as minor and major axes. We have studied only orientations in the range of 0 to 90°. However, retardation factors are symmetrical to both principal orientations due to the symmetry in network geometry and the flow domain. To indicate this symmetry the retardation ellipses have been drawn for the entire range of a (0-360°) . The good fit of these empirical retardation ellipses indicates that the influence of network orientation on retardation is systematic and predictable. No-flow boundary conditions were used for the top and bottom boundaries in our analysis to ensure that sufficient particles exit at the downgradient boundary. This procedure could potentially introduce artificial flow and transport behavior ("wall-effects"). To test for such boundary effects transport simulations were also performed for the two principal orientations in a flow domain 10x20 and 10x30, that is twice and three times the width of the flow domain shown in Figure 4.1. No significant differences were found for the various <R>ref suggesting that the observed anisotropy in retardation is not influenced significantly by artificial boundary 171 effects. The retardation ellipses serve to illustrate both the nonuniformity (different values for <R>ref) as well as the anisotropy (directional dependence of <R>ref) of retardation at the plume scale. Both aspects of retardation are governed by the variability in the aperture distribution. Anisotropy in retardation is greatest for o equal to 0 (Fig. 4.9a). A comparison of the ellipses drawn for <R>5 and <R>95 for this fracture network (CJ=0) shows that the major axis of <R>5 corresponds to the minor axis for <R>95 (Fig. 4.9a) . The degree of anisotropy (ratio of major and minor axis) for those two reference points are comparable. For high a, however, the anisotropic nature of retardation is greatly reduced (Fig. 4.9b). In other words the increased variability in fracture aperture has the effect of decreasing the directional dependence of retardation. This is true in particular for the tail of the plume (see <R>95 in Fig. 4.9b) . The opposite trends are observed for the influence of o on nonuniformity in retardation. For o equal to zero the retardation ellipses of the various reference points are relatively close to each other indicating little nonuniformity (Fig. 4.9a) . In fact, retardation is uniform for an intermediate orientation of -45°. An increase in a results in greater nonuniformity in retardation which is apparent in a clear separation of the retardation ellipses so that <R>5 is consistently lower and <R>95 is consistently greater than <R> 172 for any orientation of the network within the flow domain (Fig. 4.9b). Retardation of the mean arrival time is invariant since it is independent of the orientation of the network as well as the spread in fracture aperture. In polar coordinates, <R>m plots as a circle with its radius equal to the geometric estimate of <R> for any a (Fig. 4.9) . Note that all <R>ref of scenarios 3-4 would fall on this circle if retardation was isotropic and uniform. In theory, knowledge of the principal retardation factors for any given reference breakthrough fraction would be sufficient to predict retardation for any orientation of the network. Unfortunately, an analytical formulation of the principal retardation factors has not yet been found. It is clear, however, that the magnitude of the principal retardation factors is directly related to: (i) the aperture distribution (described by bm and (J) for each fracture set, and (ii) the network geometry (described by fracture density and fracture orientation). In the final section we address the influence of variability in fracture orientation on the extent of nonuniformity and anisotropy in plume retardation. 173 4.5.5 Influence of Variability in Fracture Orientation on Retardation In scenarios 5-6 we study retardation in networks in which the orientation of the two fracture sets (with differing mean aperture for the sets) is not strictly orthogonal. Transport was simulated in two network geometries with (i) 10° standard deviation about the two mean (orthogonal) orientations ("nearorthogonal orientation") and (ii) 90° standard deviation about the two mean (orthogonal) orientations ("random orientation") for the two fracture sets (Table 4.1). Figures 4.10a and 4.10b show example realizations of the networks with near-orthogonal and random orientation, respectively. In scenarios 5-6 migration experiments were run only in the two principal orientations (0°and 90°; see Table 4.1) for o equal to 0.0 and 0.5. We have interpolated the retardation behavior to all intermediate orientations by drawing retardation ellipses using the principal retardation factors obtained for a equal to 0.0 and 0.5 (Fig. 4.11 and 4.12, respectively). The retardation ellipses for the near-orthogonal network plot in the upper half of the diagram (Fig. 4.11a and 4.12a) and those for the random network plot in the lower half of the diagram (Fig. 4.11b and 4.12b). For the case where all apertures of one set are 80pm, and for the second set 20pm (i.e. O=0.0), the near-orthogonal network still shows directional dependence in retardation 174 (Fig. 4.11a). As might be expected, a random fracture network shows no directional dependence in retardation (Fig. 4.11b). This gradual loss in anisotropy is accompanied with an increase in the nonuniformity of retardation. The overall spread in the various <R>ref is considerably greater for the random network than for the near-orthogonal network (Fig. 4.11) . For a large spread in the aperture distribution the degree of anisotropy in retardation is greatly reduced in the nearorthogonal network (Fig. 4.12a). As a result retardation is consistently nonuniform for any orientation of the network in the flow field. Similar trends have been described already for the orthogonal network (Fig. 4.9). In the network with no preferred (i.e. random) orientation the degree of nonuniformity also increases with an increase in o (compare Fig. 4.11b and 4.12b). Based on these results the influence of variability in fracture orientation on retardation can be summarized as follows. The highest degree of anisotropy in retardation can be expected in a strictly orthogonal fracture network composed of two fracture sets with widely differing mean aperture, in particular if the spread in each aperture distribution is small (low o) . Greater variability in the orientation of the fractures for each set weakens the directional dependence of retardation. There is, however, an increase in nonuniformity, that is the overall spread in the retardation factors for 175 different reference points on the breakthrough curve. The reason for this complex dependence of retardation on the mean orientation of fracture sets in the flow field (and its variability) is again found in the coupling of fluid velocity and retardation through fracture aperture. For fracture networks with a common mean aperture (scenarios 1-2) and those fracture networks with differing mean apertures in each fracture set but no preferred mean orientation of the sets (scenario 6 ) , fluid velocity is directly proportional to b2 (i.e. the cross-sectional term in the cubic flow equation (20)). For a network in which the fracture sets have differing preferred mean orientation and aperture is nonuniformly distributed, the relationship of fluid velocity and aperture is further complicated by a correlation of local head gradients and fracture aperture. The general nature of this correlation (i.e. positive or negative) is determined by which fracture set is more closely aligned with the mean hydraulic gradient. The strength of this correlation and its influence on retardation is determined by the variability about the preferred orientation, that is, the extent to which the fracture sets are preferentially oriented. 176 4.5.6 Limitations of this Study The applicability of our findings is clearly limited by our assumptions. With respect to the fracture network we have assumed a two-dimensional network of parallel plates. In reality, a fractured medium is composed of fracture planes with rough surfaces in a three-dimensional configuration. However, the mechanisms responsible for nonuniform and anisotropic retardation are expected to be present in a threedimensional fracture network as well, namely the non-linear coupling of fluid velocity and retardation via the fracture aperture (or more precisely the specific surface area) of any given fracture segment. Direct surface profiling in natural fractures has shown that natural fractures are often rough and partially closed (e.g. Brown and Scholz, 1985). In experimental studies on flow and transport in single fractures significant channeling has been observed (e.g. Abelin et al., 1985; Brown, 1987). The parallel plate model may not be adequate to describe the channelized movement of fluid and solutes in such rough, partially-closed fractures. Moreno et al. (1988) conceptualized the fracture plane as a field of variable fracture apertures with a spatial correlation length. Using numerical migration experiments they observed channeling phenomena comparable to those observed in experimental studies. Moreno et al. (1988) also found that the retardation factor for this variable-aperture fracture could 177 be determined from the retardation equation (2) with the use of the mean aperture of the fracture. This result is consistent with our findings for retardation at the network scale. Moreno et al. did not comment on any nonuniformity in the retardation response. Based on our results we anticipate that a variability in local fracture aperture (within a single fracture plane) would also lead to nonuniform retardation even at the scale of a single fracture. At the network scale, the fracture density and the mean orientation of these rough fracture planes in three dimensions should determine to what extent this nonuniformity in retardation at the single fracture scale translates into nonuniformity and anisotropy at the fracture network scale. The general trends presented earlier regarding the influence of fracture density and variability in fracture orientation on retardation in two dimensions should in principle apply to three-dimensional networks of these rough (variable-aperture) fractures as well. Perhaps the most important assumption in the retardation approach presented here is the dependence of retardation on the surface area-to-volume ratio for a single fracture (see equation (1)). This assumption has been tested for example by Vandergraaf and co-workers who compared surface retardation factors determined from tracer BTCs in machined and natural fractures in granite to those calculated from Ka values determined from static sorption experiments using rock 178 coupons. For migration experiments in machined fractures with residence times in the order of days, retardation factors were in good agreement with those calculated from static Ka values (Vandergraaf et al., 1988). In another laboratory experiment, Vandergraaf and Drew (1991) determined a surface retardation factor for cesium in a natural fracture (0.9m by 0.9m) in a block of granite. The observed surface retardation factor (based on retardation of the 50 percent mass breakthrough fraction) was approximately half of that expected from the static Ka value. The observed differences in <R>50 and <R> could have resulted from nonuniform retardation due to the roughness of the fracture plane. Unfortunately, retardation of the mean arrival time was not reported. The surface retardation model has also been tested for transport of the strongly-sorbing uranium in constant-aperture fractures in steel (Wels, 1995). The spatial distribution of the radionuclide sorbed to the oxidized fracture walls suggested that retardation of uranium, was inversely proportional to fracture aperture, providing qualitative support for the definition of the surface retardation factor (2). These laboratory results are generally encouraging and suggest that the retardation equation (2) for a single fracture may be applicable. The retardation equation (2) for a single fracture as well as its counterpart for a fractured medium (11) assumes "ideal" sorption, that is, fast, reversible sorption that follows a 179 linear isotherm. Non-ideal sorption behavior may significantly influence retardation in single fractures and at the network scale. For example, laboratory migration experiments of strontium in constant-aperture fractures in granite showed that the influence of fracture aperture on strontium retardation was much greater than predicted by the surface retardation model (2) (Wels, 1995) . The author suggested that hysteresis in sorption, in conjunction with limited transverse mixing caused an increase in the sorption strength Ka with decreasing fracture aperture. The resulting non-linear relationship between the surface retardation factor and specific surface area would significantly increase the complexity of retardation in networks of fractures. It should be emphasized again that we have considered only surface sorption onto the rock walls in this chapter. Numerous studies have shown that retardation of a sorbing solute due to matrix diffusion and matrix sorption is potentially significant (e.g., Neretnieks et al., 1982/ Sudicky and McLaren, 1992/ Maloszewski and Zuber, 1993). Matrix sorption has a similar effect on the retardation of the solute plume as has surface sorption in most fracture networks, namely an increase in dispersion of the sorbing solute plume relative to a conservative solute plume. It may be difficult to distinguish between those two mechanism based on the comparison of conservative and reactive BTCs in a real tracer experiment. 180 4.6 Conclusions Numerical migration experiments have been performed to study retardation of reactive solutes due to surface sorption in relatively dense two-dimensional fracture networks composed of two fracture sets. The fractures are approximated by parallel plates and retardation in a single fracture is assumed to vary with the surface area-to-volume ratio of the fracture. Given these assumptions, retardation at the plume scale is a nonuniform and anisotropic process. Different segments of the plume, or equivalently different breakthrough fractions at a downstream boundary, are retarded to a different degree. The mean arrival time is the only breakthrough fraction that is consistently described by a uniform, isotropic retardation factor which is based on the specific surface area of the network. The degree to which various breakthrough fractions are retarded varies as a function of the orientation of the mean hydraulic gradient relative to the orientation of the fracture sets. This variation can be described in the form of a retardation ellipse. The degree of nonuniformity and anisotropy in the retardation factor is controlled largely by the difference in the mean apertures of each fracture set, the standard deviation in fracture aperture, and the network geometry (fracture density and fracture orientation). A larger 181 variability in fracture aperture, or in the orientation of fractures within each set, promotes nonuniform retardation while reducing the degree of anisotropy in the retardation factor. Dense fracture networks are more likely to exhibit nonuniform retardation at the plume scale than are sparse networks. The complex retardation behavior is a result of the non-linear coupling of fluid velocities and retardation via fracture aperture. Our results suggest that a transport model based on the conventional advection-dispersion equation, using a uniform retardation factor, may be conceptually incorrect, even for a dense fracture network. While the effects of nonuniform retardation modify the dispersive flux in a way that could be described by "effective" dispersion coefficients, it is our opinion that a model formulation is needed that clearly separates chemical effects from those due to hydrodynamic dispersion. One alternative approach to model the process of retardation may be the stochastic continuum approach to transport developed by Schwartz and Smith (1988) . In this approach, transport is modeled by repeatedly sampling from velocity distributions uniquely defined for different orientations of solute movements. A consideration of retardation would involve modifying these velocity distributions based on discrete network simulations such as discussed in this chapter. In this way the nonuniform and anisotropic nature of retardation at the network scale is 182 accounted for in the model. 4.7 Acknowledgements We thank Roger Beckie and Tom Clemo for the many discussions which have been helpful in formulating our ideas. This research has been supported by a grant from the Natural Science and Engineering Research Council of Canada (NSERC). 183 4.8 References Abelin, H., I. Neretnieks, S.Tunbrant, and L. Moreno, Final report of the migration in a single fracture - Experimental results and evaluation, Stripa Project IR 85-03, 1985. Bear, J., The dynamics of fluids in porous media. 7 64 pp. Elsevier, New York, 1972. Bellin, A., A. Rinaldo, W.J.P. Bosma, S.E.A.T.M. van der Zee, and Y. Rubin, Linear equilibrium adsorbing solute transport in physically and chemically heterogeneous porous formations, 1. Analytical solutions, Water Resour. Res., 29, 4019-4030, 1993. Bosma, W.J.P., A. Bellin, S.E.A.T.M. van der Zee, and A. Rinaldo, Linear equilibrium adsorbing solute transport in physically and chemically heterogeneous porous formations, 2. Numerical results, Water Resour. Res., 29, 4031-4043, 1993. Brown, S.R. and C.H. Scholz, Broad bandwidth study of the topography of natural rock surfaces, J. Gephys. Res. 90 (B14) : 12575-12582, 1985. Brown, S.R., Fluid flow through rock joints: the effect of surface roughness, J. Gephys. Res. 92(B2):1337-1347, 1987. Brusseau, M.L. and P.S.C. Rao, Sorption non-ideality during organic contaminant transport in porous media, Crit. Rev. in Env. Control 19(1), 33-99, 1989. Burkholder, H.C., Methods and data for predicting nuclide migration in geoogic media, Intern Symp. on Management of Wastes from the LWR Fuel Cycle , Denver, Colorado, 1976 Chernyshev, S.N. and W.R. Dearman, Rock Fractures, Butterworth-Heinemann, London, 1991. Clemo, T., Dual permeability modeling of fractured media, Ph.D. Thesis; University of British Columbia, 1994. Cvetkovic, V.D. and A.M. Shapiro, Mass arrival of sorptive solute in heterogeneous porous media, Water Resour. Res., 26, 2057-2067, 1990. Dverstorp, B., J. Andersson, and W. Nordqvist, Discrete fracture network interpretation of field tracer migration in sparsely fractured rock, Water Resour. Res. 28, 2327-2343, 1992. Freeze, R.A. and J.A. Cherry, Groundwater, Prentice Hall, N.J. 604 pp. 1979. 184 Garabedian, S.P.L., L. Gelhar, and M.A. Celia, Large-scale dispersive transport in aquifers: Field experiments and reactive transport theory, Rep. 315, Ralph M. Parsons Laboratory, Mass. Inst. Technol., Cambridge, 1988. Hull, C.H., J.D. Miller, and T.M. Clemo, Laboratory and simulation studies of solute transport in fracture networks, Water Resour. Res., 23, 1505-1513, 1987. Maloszewski, P. and A. Zuber, Tracer experiments in fractured rocks: Matrix diffusion and the validity of models, Water Resour. Res., 29, 2723-2735, 1993. Moreno, L., Y.W. Tsang, C.F. Tsang, F.V. Hale, and I. Neretnieks, Flow and transport in a single fracture: A stochastic model and its relation to some field observations, Water Resour. Res., 27, 2033-2048, 1988. Neretnieks, I. A note on fracture flow mechanisms in the ground, Water Resour. Res., 19(2), 363-370, 1983. Neretnieks, I., T. Eriksen, and P. Tahtinen, Tracer movement in a single fissure in granitic rock: Some experimental results and their interpretation, Water Resour. Res. 18(4), 849-858, 1982. Nordqvist, A.W., Y.W. Tsang, C.F. Tsang, B.Dverstorp, and J. Andersson, A variable aperture fracture network model for flow and transport in fractured rocks, Water Resour. Res. 28(6), 1703-1713, 1992. Raven, K.G., Hydraulic characterization of a small groundwater flow system in fractured monzonite gneiss, National Hydrol.Res.Inst., Envir. Canada, Paper No. 30, 1986. Schwartz, F.W. and L. Smith, A continuum approach for modeling mass transport in fractured media, Water Resour. Res., 24(8), 1360-1372, 1988. Smith, L., T. Clemo, and M.D. Robertson, New approaches to the simulation of field-scale solute transport in fractured rocks, Proc. 5th Canadian/American Conf. on Hydrogeology, National Groundwater A s s o c , 1990. Snow, D.T., The frequency and aperture of fractures in rock, Int. J. Rock Mech. Miner. Sci. 7, 23-40, 1970. Sudicky, E.A. and R.G. McLaren, The Laplace transform Galerkin technique for large-scale simulation of mass transport in discretely fractured porous formations, Water Resour. Res., 28, 499-514, 1992. 185 Tsang, Y.W., C.F. Tsang, I. Neretnieks, and L. Moreno, Flow and tracer transport in fractured media: A variable aperture channel model and its properties, Water Resour. Res., 24, 2049-2060, 1988 Valocchi, A.J., Spatial moment analysis of the transport of kinetically adsorbing solutes through stratified aquifers, Water Resour. Res., 25, 273-279, 1989. Vandergraaf, T.T. and D.J. Drew, Laboratory studies of radionuclide transport in fractured crystalline rock, 3rd Int. Symp. on advanced nuclear energy research, Mito City, Japan, unpublished manuscript, 1991. Vandergraaf, T.T., D.M. Grondin and D.J. Drew, Laboratory radionuclide migration experiments at a scale of one meter. Mat. Res. Soc. Symp. Proc. Vol 112, p. 159-168, 1988. Wels, C , Chapter Three: The influence of specific surface area on solute transport in fractures - An experimental analysis, In: Transport of sorbing solutes in fractured media - A numerical and experimental analysis of dispersion and retardation, Ph.D. thesis, Univ. of B.C., Vancouver, 1995. 186 Table 4.1. Summary of fracture network input parameters which were varied in the seven transport experiments. All fracture networks are comprised of two sets of fractures. The mean fracture length of both sets is 1.0m in all runs. Scenario Mean Fracture Aperture (in Turn) Fracture Orientation1 (in ') Fracture Scan Line Density2 (in fract./m) b.1 b* a, a, 80 80 0.0 90.0 6.0 6.0 80 80 0.0 90.0 3.0 3.0 0.0 90.0 6.0 6.0 6.0 6.0 6.0 6.0 80 20 >0.0 0.0 80 20 near-orthogonal 0.0+10.0 90.0±10.0 (90.0110.0 0.0±10.0) 80 20 random 0.0±90.0 90.0±90.0 (90.0±90.0 0.0±90.0) 1 orientation of each fracture set with respect to the overall hydraulic head gradient (i.e. w.r.t. the horizontal); in scenarios 5-6 two (principal) orientations of the network were studied (see text) 2 fracture densities used for generation of the network; densities actually observed for the hydraulically connected network within the flow domain are slightly lower (by no more than 3% of the input value) 187 mmmmmmmr- ] g .Iglililfg Source line ->dh/dx=0.01 ho i 5 10 Domain length (m) Figure 4.1. Flow domain, boundary conditions and the source zone used in the simulations. No-flow boundaries are imposed at the top and bottom of the model domain. 188 10 JJJ 1 1 1 1 - j , 1 i—I - = l J 8 | 1 II J 1 t 1 I ^P I I : u [j 1 M ll -=LJj== ; 1 4 Mi 1 1 , .'I ' - T - l 1 PI r 1 — 4 — i 1 I— =•-!-- J HJW n , D+frTTPl —1 _ J — ^ r>n L i' 1 |'—h fl-H-H-H ... .. _ -ffl ffl i F+rUmJ rim ST n n,Hfl HHH H— IH—• r 1 LUy | N |7U ' hHltmPBd=ffi i MHH—n * J rr nhr^^tl H~^ \ rffl—H fTTHIIl £ TT =:_= 0 1 ij-jl 1 !Y~^ r i Bf+jB f-n i ft ffl p_ _in v v —| ji w 4 > y .._—— ji IjtiJ 1 y 6 1 n _ d E3 F T -.U-- r--7R- ' i " r H i t m r n n Pit • ffnltj—Tfl~H l LI 0 i i- i 21 i uXJi 4 JI llil llil-l 611—4 " H i ^8H i i \—1 10 H o r i z o n t a l Dimension (m) Figure 4.2. Single realization of a fracture network representing the reference geometry used in scenarios 1 and 34. Only the hydraulically-connected network is shown. See Table 4.1 for fracture statistics. 189 40 "i 1 1 1 1 1 r Scenario 1 30 > :—i E-H T95 20 > M en 10^-i-i- o 0 t J l L 0 I I 0.2 "i 150 1 Tm T50 T, * r Ed 0.4 i i r i 0.6 1 r .1 1 <R> 95 CD * — * r -? 100 A - i- <R>r A K V IS-. 50 r <R> 50 <R>. 0 0 0.2 0.4 0.6 a of log 1 0 (b) Figure 4.3. Transport behavior in the dense network of scenario 1. The two fracture sets have a common aperture distribution (bm=80p.m) . Shown are the average values for conservative arrival times Tref (upper panel) and retardation factors <R>ref (lower panel) . The standard deviation of these estimates is indicated by the one-sided errorbars. The geometric estimate of <R> in scenario 1 is equal to 101 (solid line). 190 10 I I I I I I I I I I I . Non—sorbing I 1 1 1 I I 1 To = 2 d a y s 8 - a o * • l-l m CD 6 *": a^" r * — o • rH +-> IH CD > 2 —- S p a t i a l Moments: — " M l = (4.35, 5.16) _ MB= (1.62, 1.15) i 0 i i 1 i i i i i i i 2 4 6 Horizontal Dimension (m) i i i 8 i i 10 Figure 4.4a. Spatial particle distributions for a single realization of scenario 1 with o equal to 0.5 for the conservative solute plume at t0=2 days. The spatial distribution of 5000 particles was discretized in a grid of 1000 by 1000 cells. The number of particles per grid cell are indicated by the number of vertices used to draw the data point. The first (Mx) and second (M2) spatial moments (x,y) of the respective solute plumes are shown in each panel as well. 191 10 I I I I I I I I I I I I I " Sorbing i i i = T *<R> Tl x i 0 1 i 50 — 154 d a y s • • •* »M» 8 '-M»JWH « •*••• PI o • r-H w - X 6 CD - ~..JP_IJ.. ^ O .....A. * ii<* * * * - % ii i—i cd 1 |i - u •r-1 — • * . " -i : T _ — - —Spatial M o m e n t s : <D > " M t = (4.75, 5.20) _ M2= (1.9B, 1.27) . 0 u 8 i i i 1 i i i 1 i a i i i i i i i i •—i i i i : Sorbing I i 1 1 i ! " 1 1 i i I 10 i i i Tg . v<H>m * • • - *v« • xJ * - * ' ; *' * — I = 202 d a y s — 1 ,g 6 CD O 1 i 2 4 6 8 Horizontal Dimension (m) o .1—1 i- i i- p. __- * ' 4 ^ -7^4 T ,':4.ii' 4.y 4 — ^s*.. j *.,4 : — CO a u CD > 2 —Spatial M o m e n t s : '.. — " M 1 = (5.25, 5.27) _ M2= (2.23, 1.3B) . i i i 1 i i I 0 1 i i i 1 i i i 1 i i i 10 2 4 6 8 Horizontal Dimension (m) Figure 4.4b/c. Spatial particle distributions for a single realization of scenario 1 with a equal to 0.5 for the sorbing solute plume at tx= to*<R>50=154 days (upper panel) and t2=to*<R>m=2 02 days (lower panel) . 192 CO 4-1 LO U O • <0 o I a L os L ^ • C\2 LTD • 3. b CO CD CO rH o o H G -H +J cd CD - P I I I I I I £ * czj 3 ^^g co 1 CD -fj •i—i CO 3. b *m o o 1—1 ivvvyyvy;/ 00 a >> agggg 05 w CD > O o X—1 - io CiO d cv i o o LO o o o o o LO saouajjnooo jo # o o o o i i i o o CO o o CD i i i o o ^ i o CO r H G o N -H — - H m CO a co <D 0 CD -H M JH <u 0 b -P -H CD <U d 0 > H G o -H tn co H -P G (U CO " H > CO U CD CD CO H CO > -^ rH CD G CO O G M a -H o o ft • P O 4-1 M o cO CD ~ -P ftfl 0) x -P M cO G G X! -P -H o • O d 4-i d C O o CD d +J > cO CO U Xi rH a cO CD CO (0 co ^ - G u a ty i i o o cv saouajjnooo jo # 0 > O C o •P J) G co M O co -P ffi egme atio wit nel) rs_l £) a C\2 I I I I I I I I -H G G ,G • LO CO d rH (0 <tf CD CO o n +J -H g CD 3 d ) M O U 4-> M cO -p G O G -P tn <o d CD o -H M G o n Pu 4-1 CO CO —' 250 u i "i XJ i I i i r Scenario 2 T95 200 E-< > 150 i—i E^ 0 100 K CO o o y 50 % 1 0.4 "k&^Jk-- 0 0.2 i 150 r T 0.6 "i r r Tm T50 T. 1 r 100 A K V & y 0 tu / <R>50 <R>, 50 0 J 0 L 0.2 0.4 0.6 (j of l o g 1 0 ( b ) Figure 4.6. Transport behavior in the sparse network of scenario 2. The two fracture sets have a common aperture distribution (bm=80jim) . Shown are the average values for conservative arrival times Tref (upper panel) and retardation factors <R>ref (lower panel) . The standard deviation of these estimates is indicated by the one-sided errorbars. The geometric estimate of <R> in scenario 2 is equal to 101 (solid line). 194 Scenario 1 PS o .1—i w m a x . flow: Pi CD a 4.6e-03 m 3 / d cti o • r-l CD > 0 10 2 i—i i i 4 6 8 10 n Pi O m Scenario 2 m a x . flow: Pi CD £ • i-H Q 1.2e-04 m 3 / d i — i cti o • r-l CD > 0 2 4 6 8 10 Horizontal Dimension (m) Figure 4.7. Flow distribution in single realizations of the dense network of scenario 1 (top) and the sparse network of scenario 2 (bottom) for o equal to 0.5. The relative thickness of the bar indicates the amount of flow (volumetric flow rate) passing through the fracture segment. 195 QJ A V 250 200 CD A K V 150 100 50 J 0 0 J I 0.2 I L I 0.4 I I 0.6 a of log 1 0 (b) Figure 4.8. Transport behavior in the dense network of scenario 3 (top) and scenario 4 (bottom). The two fracture sets have differing aperture distributions (bm1=8 0|im and bm2=20pim) . Scenarios 3 and 4 differ only in their orientation with respect to the mean hydraulic gradient; hence, the geometric estimate of <R> is equal (<R>=161) in both scenarios (solid lines). 196 <R>, (7 = 0.0 20( / EfcSL ./' 100 • I 1 1 ,1 h .1 • u \ Vl \ 1 1 1 <R>, 1 100; : <R>, cr=0.5 /200 1/ • <R> 5 • <R>50 <& < R > m A <R>g5 fn 200 <R> > Figure 4.9. Retardation factors plotted in polar coordinates for a range of network orientations in the flow domain for a equal to 0.0 (top panel) and a equal to 0.5 (bottom panel). The data points represent the average values of <R>ref observed for the four reference mass breakthrough fractions. The various lines represent retardation ellipses for the four reference points (see text). 197 lr i I L "TW ' 1 1 ^UW PI ^Tk^ 8 o • I—I w pj Scenario 5 nearorthogonal 6 CD 4 o • 1—1 •+-> 2 3TW-T- 1i r*^£-_L] • 1 •] _L til CD > -*—^ n 0 2 f-F=^ 4 6 8 10 10 Pi O 8 • i—i w Pi Scenario 6 random 6 CD B •i—i Q 4 i—i ctf o •i—i 2 CD > 0 0 2 4 6 8 10 Horizontal Dimension (m) Figure 4.10. Example realization of a near-orthogonal fracture network in scenario 5 (top panel) and a random network in scenario 6 (bottom panel). Only the hydraulicallyconnected network is shown. 198 <R>, near-orthogonal cr=0.0 random CJ=0.0 Figure 4.11. Retardation ellipses for the near-orthogonal network of scenario 5 (top panel) and the random network of scenario 6 (bottom panel) for o equal to 0.0. The ellipses are constructed based on transport simulations with the networks oriented in the two principal orientations (0° and 90°) . Retardation ellipses are drawn for the 5%, 50%, mean, and 95% breakthrough fractions. 199 <R>, near—orthogonal (7=0.5 random cr=0.5 Figure 4.12. Retardation ellipses for the near-orthogonal network of scenario 5 (top panel) and the random network of scenario 6 (bottom panel) for a equal to 0.5. The ellipses are constructed based on transport simulations with the networks oriented in the two principal orientations (0° and 90°) . Retardation ellipses are drawn for the 5%, 50%, mean, and 95% breakthrough fractions. 200 Chapter Five: CONCLUSIONS AND RECOMMENDATIONS 5.1 Concluding remarks In this thesis we have studied the transport of sorbing solutes in fractured media using numerical migration experiments in parallel-plate fractures (Chapter Two) and network of fractures (Chapter Four), as well as laboratory migration experiments in constant-aperture fractures (Chapter Three). The numerical analyses suggest that (1) surface sorption may enhance dispersion at the scale of a single fracture, and (2) retardation may be nonuniform and anisotropic at the network scale. Four key assumptions have been made in order to arrive at these conclusions: (i) a fracture can be represented by a pair of parallel plates, (ii) only solute in the immediate vicinity to the fracture wall can participate in surface sorption, (iii) sorption is a fast, reversible process which can be described by a linear isotherm, and (iv) retardation varies in proportion to the surface area-to-volume ratio of a fracture. We emphasize that the applicability of our numerical results is limited by these assumptions. The laboratory migration experiments described in 201 Chapter Three clearly illustrate this point. The results of these experiments indicate that hysteresis in sorption (violation of assumption (iii)) results in a much stronger dependence of retardation on fracture aperture than is postulated in assumption (iv). One of the primary goals of research on reactive solute transport in fractured media is to describe and predict contaminant transport in a natural environment, for example in a fractured rock formation hosting a waste site. Natural systems are more complex than has been assumed in the numerical analyses or has been engineered in the laboratory experiments. Future research is needed to show to what extent our conclusions apply to transport in complex, natural systems. It is recommended that research be carried out to investigate the sensitivity of our findings to the simplifying assumptions made in this study. To this end several research topics are identified below that should be addressed in order to improve our understanding of reactive solute transport in fractured media. In Chapter Two we assumed that only solute in close proximity to the fracture wall can participate in the sorption process. It has been shown that the localized nature of surface sorption, in conjunction with limited transverse mixing, results in enhanced longitudinal spreading of the sorbing solute relative to that of a conservative solute. These findings (and the underlying assumptions) could be 202 tested in controlled laboratory migration experiments in constant-aperture fractures. According to our results the degree of enhanced dispersion should decrease significantly in the transport regime of mixed dispersion. This hypothesis can be tested by repeating migration experiment at different fluid velocities. It is important,, however, to choose a combination of sorbing substrate and solute that exhibits fast, reversible sorption to avoid enhanced dispersion due to non-ideal sorption effects (see Chapter Three). One possibility would be to coat an impermeable fracture surface such as plexiglass with an exchange resin so that the sorption process is fast and well-defined. All transport models used in this thesis assumed a fast, reversible sorption which follows a linear isotherm. This assumption is the basis for the surface retardation model which states that retardation is proportional to the specific surface area of a fracture. Our experimental results (Chapter Three) indicate that the dependence of retardation on specific surface area may be significantly greater than proposed by this model for a solute that exhibits hysteresis in sorption. Likewise, other non-ideal sorption behaviour such as non-equilibrium or non-linear sorption can be expected to influence dispersion and retardation of a sorbing solute. The influence of hysteresis in the sorption reaction (and other non-ideal sorption effects) on retardation should be studied in more detail using numerical and laboratory migration 203 experiments. Our numerical migration experiments in networks of fractures (Chapter Four) indicate that, at the plume scale, retardation may be nonuniform and anisotropic. This retardation behaviour is a result of the cross correlation between the retardation, factor and fluid velocity in each fracture, which follows directly from the assumed dependence of the surface retardation factor and fluid velocity on fracture aperture. These assumptions appear justified in a network of smooth, open fractures which can be represented by parallel plates. It is uncertain, however, whether they apply in a network of rough-walled, partially-closed fractures. We anticipate difficulties in testing the dependence of the surface retardation factor on fracture aperture (or specific surface area) for rough-walled fractures for two reasons. First, it is difficult to measure fracture aperture (or specific surface area) accurately in-situ with existing measuring techniques. Second, the measured fracture aperture (specific surface area) of a rough-walled fracture, presumably some "average" value of the fracture plane, may not coincide with the "effective" fracture aperture "sampled" by the solute plume as a result of channeling. It may be advantageous to relate retardation to fracture parameters which are easier to measure in the field than fracture aperture (or specific surface area). Of particular interest would be the relationship between retardation and 204 transmissivity of a fracture. The surface retardation model studied in this thesis would suggest a negative correlation between the surface retardation factor and transmissivity of a fracture, assuming a wide fracture has a greater transmissivity than a narrow fracture. This relationship could be tested by measuring the retardation of a sorbing solute during transport in rough-walled fractures which differ in transmissivity. These experiments could be carried out in the laboratory or at the field scale. A negative cross correlation between retardation and transmissivity would result in nonuniform and anisotropic retardation at the network scale in much the same way as has been illustrated in Chapter Four. Direct experimental verification of this proposed retardation behaviour at the network scale may become feasible with the current development of large-scale field sites in fractured rock. It has been concluded that the use of a transport model based on the conventional advection-dispersion equation, using a uniform retardation factor, would be conceptually incorrect in a fractured medium that exhibits nonuniform and anisotropic retardation. Research should focus on the theoretical development of continuum models that accommodate nonuniform and anisotropic retardation. Such a transport model would improve our ability to describe and predict the migration of sorbing contaminants in fractured media. 205 5.2 Summary of Conclusions In Chapter Two the transport of a sorbing solute in a parallel-plate fracture was simulated using a random walk model which assumes that sorption occurs only in the immediate vicinity of the fracture wall. The numerical migration experiments suggest: (1) Solute transport in a parallel plate fracture is described accurately by the one-dimensional form of the advectiondispersion equation provided the condition of transverse homogenization is met. A sorbing solute requires a greater transport distance for transverse homogenization to be achieved than a non-sorbing solute. (2) A sorbing solute experiences greater longitudinal spreading during transport in a parallel-plate fracture than does a conservative solute. The magnitude of this enhanced dispersion increases with higher fluid velocities and reaches a maximum in the range of fluid velocities characteristic of Taylor dispersion. An effective dispersion coefficient has been developed which accounts for the additional longitudinal spreading of a sorbing solute plume in a parallel-plate fracture. 206 In Chapter Three, the sorbing radionuclides strontium and uranium were injected in smooth-walled, constant-aperture fractures in granite and steel, respectively. The laboratory migration experiments indicate: (3) The retardation of uranium was inversely proportional to fracture aperture providing qualitative support for the definition of the surface retardation factor. (4) The influence of fracture aperture on the retardation of strontium was much greater than predicted by the surface retardation factor. Strontium retardation was approximately an order of magnitude greater in a narrow fracture (b=450}im, Ra=45) compared to in a wide fracture (b=780jim, Ra=3.5) . It is hypothesized that hysteresis in sorption, in conjunction with limited transverse mixing across the aperture, caused the apparent increase in sorption strength with a decrease in fracture aperture. (5) The reactive solute strontium experienced significantly more dispersion than the non-reactive tracer tritium. Model simulations suggested that the dispersion of strontium was approximately an order of magnitude greater than that due to hydrodynamic dispersion. This enhanced dispersion is believed to be a result of chemical heterogeneity, hysteresis in sorption, and limited transverse mixing. 207 In Chapter Four, the movement of a sorbing solute at the scale of a fracture network was simulated assuming the surface retardation model holds at the scale of a single fracture. The numerical migration experiments suggest: (6) Retardation at the plume scale is a nonuniform and anisotropic process. Different segments of the plume, or equivalently different breakthrough fractions at a downstream boundary, are retarded to a different degree. The degree to which various breakthrough fractions are retarded varies as a function of the orientation of the mean hydraulic gradient relative to the orientation of the fracture sets. This variation can be described in the form of a retardation ellipse. (7) The degree of nonuniformity and anisotropy in the retardation factor is controlled largely by the difference in the mean apertures of each fracture set, the standard deviation in fracture aperture, and the network geometry (fracture density and fracture orientation). A larger variability in fracture aperture, or in the orientation of fractures within each set, promotes nonuniform retardation while reducing the degree of anisotropy in the retardation factor. 208 (8) The mean arrival time is the only breakthrough fraction that is consistently described by a uniform, isotropic retardation factor which is based on the specific surface area of the network. It is concluded that a transport model based on the conventional advection-dispersion equation, using a uniform retardation factor, may be conceptually incorrect, even for a densely fractured medium. The present analysis of reactive solute transport in fractures and network of fractures suggests that the concept of retardation, originally developed for porous media, may require modifications before it can be applied to a fractured medium. To this end future research should be conducted to (i) test experimentally the validity of the surface retardation model in single fractures or channels for solutes and transport conditions of interest, (ii) study the influence of non-ideal sorption behaviour on retardation and dispersion of solutes in fractured media, (iii) examine the relationship of retardation and transmissivity in rough-walled fractures, and (iv) develop a continuum transport model that accommodates nonuniform and anisotropic retardation. 209
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Transport of sorbing solutes in fractured media : a...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Transport of sorbing solutes in fractured media : a numerical and experimental analysis of dispersion… Wels, Christoph 1995
pdf
Page Metadata
Item Metadata
Title | Transport of sorbing solutes in fractured media : a numerical and experimental analysis of dispersion and retardation |
Creator |
Wels, Christoph |
Date Issued | 1995 |
Description | The primary pathways for contaminant transport in lowpermeability fractured rock are likely to be through a network of hydraulically-connected fractures in the rock formation. Sorption of contaminants to the fracture walls may significantly retard their transport. The influence of sorption on solute transport is analyzed using both numerical and laboratory migration experiments. A random walk model is developed to simulate solute transport in a parallel plate fracture, assuming that fast, reversible, and linear sorption occurs in a small zone adjacent to the fracture wall. It is shown that a sorbing solute experiences greater longitudinal spreading than does a conservative solute. The magnitude of this enhanced dispersion reaches a maximum in the range of fluid velocities characteristic of Taylor dispersion. At the network scale, transport is simulated by tracking particles in discrete fracture networks, assuming that within each fracture, retardation varies in proportion to the product of a surface distribution coefficient and the specific surface area of the fracture. The results suggest that retardation at the plume scale is a non-uniform and anisotropic process. Different segments of the plume, or equivalently different breakthrough fractions at a downstream boundary, are retarded to a different degree. The degree to which various breakthrough fractions are retarded varies as a function of the orientation of the mean hydraulic gradient relative to the orientation of the fracture sets. This variation can be described in the form of a retardation ellipse. To test the surface retardation model used in the numerical analyses, the sorbing radionuclides strontium and uranium were injected in smooth-walled fractures in granite and steel, respectively. The retardation of uranium was inversely proportional to fracture aperture, providing qualitative support for the definition of the surface retardation factor. The influence of fracture aperture on the retardation of strontium was much greater than that predicted by the surface retardation factor. Strontium retardation was approximately an order of magnitude greater in a narrow fracture (b=450jim, Ra=45) compared to that in a wide fracture (b=780jim, Ra=3.5) . It is hypothesized that hysteresis in sorption, in conjunction with limited transverse mixing across the fracture, caused the apparent increase in sorption strength with a decrease in fracture aperture. |
Extent | 12374911 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052914 |
URI | http://hdl.handle.net/2429/7459 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geological Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1995-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1995-060810.pdf [ 11.8MB ]
- Metadata
- JSON: 831-1.0052914.json
- JSON-LD: 831-1.0052914-ld.json
- RDF/XML (Pretty): 831-1.0052914-rdf.xml
- RDF/JSON: 831-1.0052914-rdf.json
- Turtle: 831-1.0052914-turtle.txt
- N-Triples: 831-1.0052914-rdf-ntriples.txt
- Original Record: 831-1.0052914-source.json
- Full Text
- 831-1.0052914-fulltext.txt
- Citation
- 831-1.0052914.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052914/manifest