SOLUTE TRANSPORT THROUGH FRACTURED ROCKS - DiVA

SOLUTE TRANSPORT THROUGH FRACTURED ROCKS: THE INFLUENCE OF GEOLOGICAL HETEROGENEITIES AND STAGNANT WATER ZONES

Batoul Mahmoudzadeh Doctoral Thesis in Chemical Engineering

Department of Chemical Engineering and Technology School of Chemical Science and Engineering Royal Institute of Technology Stockholm, Sweden 2016

Solute transport through fractured rocks: the influence of geological heterogeneities and stagnant water zones Batoul Mahmoudzadeh Doctoral Thesis in Chemical Engineering TRITA-CHE Report 2016:15 ISSN 1654-1081 ISBN 978-91-7595-906-1

KTH School of Chemical Science and Engineering SE-100 44 Stockholm SWEDEN

Akademisk avhandling som med tillstånd av Kungl Tekniska Högskolan framlägges till offentlig granskning för avläggande av teknologie doktorsexamen i kemiteknik tisdagen den 26 april 2016 klockan 10.00 i Sal F3, Lindstedtsvägen 26, Kungliga Tekniska Högskolan, Stockholm. Avhandlingen försvaras på engelska. Fakultetsopponent: Dr. David Hodgkinson, Ordförande i Quintessa.  Batoul Mahmoudzadeh 2016 Tryck: US-AB

ii

To my parents

iii

iv

Abstract To describe reactive solute transport and retardation through fractured rocks, three models are developed in the study with different focuses on the physical processes involved and different simplifications of the basic building block of the heterogeneous rock domain. The first model evaluates the effects of the heterogeneity of the rock matrix and the stagnant water zones in part of the fracture plane. It accounts for the fact that solutes carried by the flowing water not only can diffuse directly from the flow channel into the adjacent rock matrix composed of different geological layers but can also at first diffuse into the stagnant water zones and then from there further into the rock matrix adjacent to it. The second and the third models are dedicated to different simplifications of the flow channels. Both account for radioactive decay chain, but consider either a rectangular channel with linear matrix diffusion or a cylindrical channel with radial matrix diffusion. Not only an arbitrary-length decay chain, but also as many rock matrix layers with different geological properties as observed in the field experiments can be handled. The analytical solutions thus obtained from these three models for the Laplace-transformed concentration in the flow channel can all be conveniently transformed back to the time domain by use of e.g. de Hoog algorithm. This allows one to readily include the results into a fracture network model or a channel network model to predict nuclide transport through flow channels in heterogeneous fractured media consisting of an arbitrary number of rock units with piecewise constant properties. More importantly, the relative impacts and contributions of different processes in retarding solute transport in fractured rocks can easily be evaluated by simulating several cases of varying complexity. The results indicate that, in addition to the intact wall rock adjacent to the flow channel, the stagnant water zone and the rock matrix adjacent to it may also lead to a considerable increase in retardation of solute in cases with narrow channels. Moreover, it is elaborated that, in addition to the chain decay, it is also necessary to account for the rock matrix comprising two or more geological layers in safety and performance assessment of the deep repositories for spent nuclear fuel. The comparison of the simulation results with respect to either rectangular or cylindrical channels further suggest that the impact of channel geometry on breakthrough curves increases markedly as the transport distance along the flow channel and away into the rock matrix increase. The effect of geometry is more pronounced for transport of a decay chain when the rock matrix consists of a porous altered layer. Additionally, a model is developed to study the evolution of fracture aperture in crystalline rocks mediated by pressure dissolution and precipitation. It accounts for not only advective flow that can carry in or away dissolved minerals but also the fact that dissolved minerals in the fracture plane, in both the flow channel and the stagnant water zone, can diffuse into the adjacent porous rocks. The analytical solution obtained in the Laplace space is then used to evaluate the evolution of the fracture aperture under combined influence of stress and flow, in a pseudo-steady-state procedure. The simulation results give insights into the most important processes and mechanisms that dominate the fracture closure or opening under different circumstances. It is found that the times involved for any changes in fracture aperture are very much larger than the times needed for concentrations of dissolved minerals to reach steady state in the rock matrix, the stagnant water zone and the flow channel. Moreover, it is shown that diffusion into the rock matrix, which acts as a strong sink or source for dissolved minerals, clearly dominates the rate of concentration change and consequently the rate of evolution of the fracture.

v

Sammanfattning Tre modeller har utvecklats för att beskriva hur radionuklider och andra lösta ämnen transporteras och fördröjas under inverkan av olika fysikaliska processer när de förs av sipprande vatten i kanaler i sprickor i berg. De tre modellerna lägger tyngden på olika processer och utgör byggstenar i mer omfattande modeller där kanalerna sammanlänkas till komplicerade nätverk i heterogent berg. Den första modellen behandlar hur heterogeniteter i bergmatrisen och närvaron av stagnant vatten i sprickan vid sidan av vatten i den strömmande kanalen inverkar på transporten. Närvaron av stagnanta vattenområden kan påverka transporten genom att ämnena kan diffundera in i det stillastående vattnet och därigenom fördröjas. Dessutom kan ämnena diffundera vidare in i den porösa bergmatrisen inte bara från det strömmande vatten utan även från det stillastående och fördröjas ytterligare. Den andra och tredje modellen behandlar förenklade specialfall av kanalgeometrin, den ena med rektangulärt den andra med cirkulärt tvärsnitt. I båda fallen behandlas radionuklider som har olika egenskaper och hela sönderfallskedjor. Modellerna kan hantera kanaler med många lager av sekundära mineral som kan ha bildats på sprickvägen i såväl kanalen som i den stagnanta området. De analytiska lösningarna till modellerna i har erhållits med hjälp av Laplacetransformer och återtransformeras numeriskt till tidsrymden med hjälp av de Hoog–algoritmen. Detta gör det möjligt att enkelt bygga upp omfattande nätverk av kanaler med olika egenskaper i komplexa tredimensionella kanalnätverk för att simulera hur kedjor av radionuklider transporteras med vatten i komplicerade heterogena bergvolymer. Varje kanal eller del av kanal ges konstanta egenskaper men kan förfinas godtyckligt i kortare kanaldelar om detta skulle krävas för att beskriva kanaler med varierande egenskaper längs kanalen. Genom att analytiska lösningar till modellerna tagits fram kan man enkelt se sambandet mellan olika processer och se hur variationer i parametrar kan påverka transporten. Ett antal simulerade fall tyder på att upptag och sorption i den porösa bergmatrisen, både från kanaler med strömmande vatten och från det stagnanta vatten har mycket stor inverkan på fördröjning av radionukliderna. Inverkan av stagnanta vattnet är särskilt stor då bredden på kanalen med det strömmande vattnet är liten. Simuleringsresultaten tyder även på att egenskaperna hos och tjocklekarna av de sekundära mineralen kan ha stor inverkan på fördröjningen av viktiga radionuklider. Jämförelser mellan cirkulära och rektangulära kanaler visar att den första geometrin visar större fördröjning än den andra ju längre transportvägen och därmed inträngningen i bergmatrisen blir vid i övrigt lika förhållanden. Närvaron av en mer porös zon på sprickväggen har olika inverkan på nuklidkedjors fördröjning i kanaler med olika form på strömningstvärsnittet. En ytterligare studie presenters där upplösning och återutfällning av de mineral som bygger upp bergmatrisen kan påverka sprickvidden och hur detta påverkas av de bergspänningar som sprickorna utsätts för. Detta är av intresse för de långa tider, hundratusentals år, som måste beaktas vid säkerhetsanalys av ett slutförvar för utbränt kärnbränsle. Samma transportprocesser och mekanismer som i de tre tidigare modellerna beaktas och samma matematiska teknik används för att lösa ekvationerna. Simuleringarna visar att det tar mycket längre tider än de nyss nämnda för att dessa processer nämnvärt skall påverka sprickaperturen och därmed nuklidtransporten. vi

Acknowledgements This work has been carried out at the Division of Chemical Engineering, Department of Chemical Engineering and Technology, of the Royal Institute of Technology, KTH, in Stockholm, Sweden. I wish to sincerely thank my colleagues and the staff for all their help and support during these years. In particular, I would like to express my deepest appreciation to the supervisor of my research project, Associate Professor Longcheng Liu, for sharing with me his great knowledge and scientific experience, and for his continuous commitment, support, guidance and encouragement not only in research work but also in many aspects of my life from the first day I came to Stockholm. Thank you very much for all your efforts to make this step a reality. Special thanks to my supervisor, Professor Ivars Neretnieks, for so many interesting and enlightening discussions throughout the research work. Your knowledge and encouragement always gave me driving force not only for the research but also for the life. I highly appreciate the good times and the interesting discussions in the social events. My sincere gratitude also goes to Professor Luis Moreno, for his much valuable advices during these years, scientifically and non-scientifically. I appreciate very much the time you have spent for questions, discussions and solutions to the problems. The Swedish Nuclear Fuel and Waste Management Company (SKB) is gratefully acknowledged for the financial support to this project. It has been very interesting and stimulating to participate in this project, and for that I would like to thank especially Björn Gylling. Thanks to all the staff and colleagues at the divisions of Chemical Engineering and Transport Phenomena. Many thanks to Associate Professor Joaquin Martinez and Associate Professor Kerstin Forsberg for all their help and consideration. Big thanks to Ann Ekqvist, Helene Hedin, Vera Jovanovic for always being friendly to give kind assistance in all administrative matters. My last, but not least, words belong to my dearest ones, my wonderful parents, for their endless love, support and encouragement throughout my life. I would like to express my deepest gratitude for all your understanding and patience. I love you and dedicate this thesis to you.

vii

List of papers included in the thesis I. Solute transport in fractured rocks with stagnant water zone and rock matrix composed of different geological layers–Model development and simulations Batoul Mahmoudzadeh, Longcheng Liu, Luis Moreno and Ivars Neretnieks Water Resources Research, 49: 1709–1727 (2013). Contribution: Main author, formulating the model and finding the analytical solutions, making the simulations and discussing the results.

II. Solute transport in a single fracture involving an arbitrary length decay chain with rock matrix comprising different geological layers Batoul Mahmoudzadeh, Longcheng Liu, Luis Moreno and Ivars Neretnieks Journal of Contaminant Hydrology, 164: 59–71(2014). Contribution: Main author, formulating the model and finding the analytical solutions, making the simulations and discussing the results.

III. Solute transport through fractured rock: radial diffusion into the rock matrix with several geological layers for an arbitrary length decay chain Batoul Mahmoudzadeh, Longcheng Liu, Luis Moreno and Ivars Neretnieks Journal of Hydrology, 536: 133–146 (2016). Contribution: Main author, formulating the model and finding the analytical solutions, making the simulations and discussing the results.

IV. Evolution of fracture aperture mediated by dissolution in a coupled flow channel– rock matrix–stagnant zone system Batoul Mahmoudzadeh, Longcheng Liu, Luis Moreno and Ivars Neretnieks Submitted for publication (2015). Contribution: Main author, formulating the model and finding the analytical solutions, making the simulations and discussing the results.

These papers are appended at the end of the thesis.

Paper I: © 2013 American Geophysical Union. Paper II: © 2014 Elsevier B.V. Paper III: © 2016 Elsevier B.V.

viii

Table of contents Abstract ...................................................................................................................................... v Sammanfattning ........................................................................................................................ vi Acknowledgements ................................................................................................................... vii List of papers included in the thesis ........................................................................................ viii 1 2

3

4

5

Introduction....................................................................................................................... 1 Transport of a stable nuclide ........................................................................................... 7 2.1 Mathematical model................................................................................................. 8 2.2 Solution in the Laplace domain ............................................................................. 13 2.3 Simulations and discussion .................................................................................... 15 Transport of nuclides in an arbitrary length decay chain .......................................... 18 3.1 Mathematical model............................................................................................... 19 3.2 Solution in the Laplace domain ............................................................................. 21 3.3 Simulations and discussion .................................................................................... 23 Evolution of Fracture Aperture..................................................................................... 28 4.1 Mathematical model............................................................................................... 30 4.2 Fracture closure rate ............................................................................................... 34 4.3 Simulations and discussion .................................................................................... 35 Concluding Remarks ...................................................................................................... 39

Notations ................................................................................................................................. 42 References ............................................................................................................................... 47 Appended papers

ix

As far as the laws of mathematics refer to reality, they are not certain, and as far as they are certain, they do not refer to reality. Albert Einstein

x

1 INTRODUCTION

In many countries, high-level radioactive wastes will be deposited in deep geological repositories in crystalline rocks (Li and Chiou, 1993; Neretnieks, 1993; Gylling et al, 1998). Since the canisters encapsulating the waste may eventually leak due to defect, corrosion or any other reasons, it is vital to develop models that may be used to efficiently describe the process of reactive transport of radionuclides through fractured rocks. The aim is to aid in safety and performance assessment of repositories for spent nuclear fuel and other radioactive wastes. Consequently, many models have been developed over recent years that simplify the fracturematrix system to different extents but all account for the most important mechanisms governing nuclide transport in fractured rocks. The differences between these modeling approaches have been discussed by Selroos et al. (2002). It was emphasized, in particular, that the discrete fracture network model, DFN, is based essentially on the assumption that some fractures are fully open over their entire size, and that they form a conductive network by intersecting each other (Hartley and Roberts, 2012; Dershowitz et al., 1999). However, the premise that rock fractures are either fully open or closed is not supported by field observations. Instead, these observations suggested that fracture surfaces are uneven, and groundwater flows mainly through open parts of fractures whereas there is a negligible flow through the rock matrix (Birgersson et al., 1992, 1993; Abelin et al., 1991, 1994; Tsang and Neretnieks, 1998; Stober and Bucher, 2006). The narrow conductive parts will form “flow channels” (Tsang and Neretnieks, 1998) and subsequently a network for water flow and solute transport from the repository to the biosphere (Neretnieks, 1993). Based on this concept, the channel network model, CNM, has been developed and used for solute transport simulations in fractured rocks (Gylling, 1997; Gylling et al., 1999; Moreno and Neretnieks, 1993a; Moreno et al., 1997). It accounted for the fact that the low permeability porous rock matrix plays a dominant role on retarding solute transport. The solutes, which are carried by the flowing water in the channel, can diffuse in and out of the immobile groundwater in the rock matrix, through “matrix diffusion” process that makes a much larger water volume accessible for the nuclides to reside in (Neretnieks, 1980; Sudicky and Frind, 1984). As a result, sorbing nuclides would find more mineral surfaces in the matrix to be attached and many nuclides could be retarded sufficiently to decay to insignificant concentrations. Both DFN and CNM models and many other models (Poteri and Laitinen, 1999; Hartley, 1998; Dershowitz et al., 1998) did not, however, consider that in heterogeneous fractured media where the fluid velocity field varies, there are some regions in the fracture plane, in which 1

advective transport is essentially negligible. The formation of these low velocity regions may be because the fracture has a small aperture there or because it contains more particle fragments forming a porous medium with low permeability. We call such regions with much less flow, “stagnant water zones”. Wherever present, the solutes carried by the flowing water in the channel may first diffuse into the stagnant water zones in the fracture plane and then from there diffuse into the porous rock matrix. For this reason, Neretnieks (2006) addressed the influence of the stagnant water zone on solute transport in fractured rocks. It was found that diffusion into the stagnant water zone could significantly retard solute transport, especially when the zone is wide. This gives the solute access to a larger surface over which diffusion into the rock matrix takes place. As a result, considerable amounts of nuclides may be retained for a long time in the rock matrix adjacent to the stagnant water zone, instead of directly diffusing into the rock matrix adjacent to the flow channel. The diffusion into the stagnant water zones in the fracture plane can be accounted for in the CNM, but not in the DFN models because of its main assumption of fully open or closed fractures. Furthermore, in most of the models for flow and solute transport through fractured rocks, it is assumed that the rock matrix is a homogenous medium with either infinite or finite thickness (Neretnieks, 1980; Moreno and Neretnieks, 1993b; Moreno et al., 1996; Barten, 1996; Grisak and Pickens, 1980; Tang et al., 1981; Sudicky and Frind, 1982; Davis and Johnston, 1984; van Genuchten, 1985; Zimmerman and Bodvarsson, 1995; Sharifi Haddad et al., 2013). However, the structure of the rock matrix, in reality, is more complex as the rock may be altered near the fracture (Löfgren and Sidborn, 2010b; Cvetkovic, 2010; Piqué et al., 2010). As a result, it is anticipated that the altered parts of the rock may strongly affect the transport of solutes in fractured media because they have significantly different retardation capacities than the intact wall rock (Selnert et al., 2009). The geological materials in the altered zone could, as shown in Figure 1, be made up of fault gouge, fracture coating, fault breccia, mylonite, cataclasite, and altered rock etc. (SKB, 2003). The altered rock appears mostly to be the result of brittle reactivation of old ductile/semi-ductile deformation zones, while the fault breccia and fault gouge originate mainly from the movements along fault planes and therefore are distributed in variable amounts and proportions over the fracture and structure planes (Löfgren and Sidborn, 2010a; Löfgren and Sidborn, 2010b). In spite of this knowledge, however, it is difficult to identify all above-mentioned matrix layers with distinctively different geological properties in field observations. For this reason, it is common to simplify the structure of the altered zone somehow and consider it as a mixture of some layers. In line with this approach, Tsang and Doughty (2003) conceptualized the rock matrix as an assemblage of the intact wall rock and several layers containing the altered material with an enhanced porosity, in order to account for their differences in geological properties. Moreno and Crawford (2009) also considered the rock matrix as being composed of several geological layers in modeling solute transport in fractured rocks. They found that the layers close to the fracture surface might be important in determining the behavior of tracer retardation during site characterization. The more deeply located layers in the rock matrix may, however, have a large impact on the nuclide transport under the conditions prevailing at the timescales of the performance assessment.

2

Figure 1. A sketch of a typical conductive structure of different rock layers, based on conceptual model taken by SKB (2002a).

With this understanding, the work presented in the thesis attempts to get an insight into the synchronous and the overall effects of both the layered rock matrix and the stagnant water zone in retarding solute transport. The constructed model combines, in essence, the two models developed by Moreno and Crawford (2009) and by Neretnieks (2006), respectively, based on an idealization of the basic building block of a heterogeneous rock domain. It accounts for not only the stagnant water zone with a finite width lying in the fracture plane but also the rock matrix, as schematically shown in Figure 2, composed of both the intact wall rock with a finite thickness and the altered zone that may be made up of altered rock, cataclasite, fault gouge and fracture coating etc. all at the same time. In addition, the model considers that the rock matrices adjacent to both the flow channel and the stagnant water zone may have different structures and geological properties. The effects of physical and chemical heterogeneities, as studied by Deng et al. (2010), can thus be taken into account by considering the fractured rocks as a series of blocks with piece-wise constant geological properties.

3

Fracture coating I

Flowpath Fracture coating II Fault gauge II

Fault gauge I

Cataclasite II

Cataclasite I

Altered rock II

Altered rock I

Intact wall rock II

Intact wall rock I

Figure 2. Schematic conceptualization of two layered rock volumes.

In addition, it is noticed that the diffusion in the rock matrix is assumed, in most models, to be one dimensional and orthogonal to the fracture surface where the fracture is conceptualized as a rectangular channel. This is true when a wide channel is involved in the system. However, since the flow channels in fractures are most possibly narrow, the matrix diffusion should be more realistically modeled as radial diffusion from a cylindrical channel, especially when the diffusion penetration depth grows larger than the width of the channel. Moreover, it should be stressed that, although fully cylindrical channels are seldom found in nature, there are some cases where cylinder-like channels may be more adequately used for modeling fluid flow and solute transport. For instance, circular layered rock formations might be built up during the completion (rock cracking, mud contamination etc.) of a horizontal well used to inject the waste. Such cylinder-like channels can also be found at fracture intersections in fractured rocks. Around these channels it may be expected that over time, chemical changes will lead to development of circular matrix layers when the diffusion/reaction penetration depth becomes larger than the channels width. As a consequence, some theoretical and experimental studies have been performed to describe the radial diffusion of solutes into a soil or rock matrix from a cylindrical macropore or fracture for different boundary conditions (van Genuchten et al., 1984; Rasmuson and Neretnieks, 1986; Carrera et al., 1998; Rahman et al., 2004; Cihan and Tyner, 2011; Chopra et al., 2015). The radial diffusion has been shown to be more effective than linear diffusion in retarding solute transport through a channel with the same flow-wetted surface (Rasmuson and Neretnieks, 1986). It has also been shown that radial diffusion gives access to a larger rock matrix volume than linear diffusion and thus causes larger retardation (Neretnieks, 2015). The developed model is then extended to study reactive transport through either a rectangular or a cylindrical channel with the rock matrix composed of several geological layers with different properties. The aim is to explore the relative impact and contribution of the altered rock matrix and the decay chain on transport of nuclides through a rectangular/cylindrical channel with linear/radial matrix diffusion. For simplicity, however, the model ignored the presence of the stagnant water zones in the fracture plane at this stage. Nevertheless, the analytical solutions, obtained in the Laplace domain, are applicable to the cases involving an arbitrary-length decay chain, where the rock matrix consists of any number of different geological layers. These analytical solutions can be used for prediction of nuclide transport through complex networks of channels, and serve as verification tools for the numerical codes. With the use of the developed models, the simulation results obtained from the rectangular 4

channel with linear matrix diffusion are also compared with those obtained from the cylindrical channel with radial matrix diffusion. As mentioned above, the flow and transport in a fractured rock mass was conceptualized in the CNM to take place in a three-dimensional network of channels; the solute from one channel enters one or more channels at an intersection. This channel network model obviates the need to account explicitly for “hydrodynamic dispersion” as the conventional advection-dispersion model does. Instead, as discussed by Neretnieks (1983), hydrodynamic dispersion may be more adequately considered to be the result of different residence times in different pathways. This eliminates the problem of the inconsistency between the use of the advection-dispersion equation and the results obtained from field experiments; the advection-dispersion equation requires a constant value of the “dispersion coefficient” whereas the observed “dispersion coefficient” was shown to increase with distance in field experiments. In addition, it is also found in numerous simulations with the channel network model that the large scale dispersion in a network is dominated by velocity dispersion and matrix diffusion; and that conventional “hydrodynamic dispersion” in individual channels has a small impact (Moreno and Neretnieks, 1991, 1993b; Gylling, 1997). For this reason and also because hydrodynamic dispersion in channels in fractures is not well understood (Molz, 2015; Tsang et al., 2015), this mechanism was not considered in the models developed for solute transport through individual channels in the study.

Pressure dissolution

Free-face dissolution/Precipitation Flow

Figure 3. Schematic view of main processes studied in this paper in control of fracture aperture evolution (left); top view of fracture plane, showing preferential flow path and stagnant water zones wherein fine-grained crystals are located (right).

In the models, we commonly assumed that the fracture is made of two parallel surfaces, meaning that the fracture aperture is constant all the time. However, rock fractures may close or open due to mechanical and chemical effects, which may considerably influence water flow and solute transport in fractured rocks (Jing et al., 2013; Zhao et al., 2014). In particular, dissolution and precipitation of minerals on fracture walls may either seal or widen a fracture. Evolution of the fracture aperture over time in stressed natural fractures has been experimentally observed (Polak et al., 2003; Beeler and Hickman, 2004, Liu et al., 2006; Yasuhara et al., 2006a), and the main causes are still under discussion. Evolution of fracture permeability has also been studied theoretically and modeled by accounting for dissolution and 5

precipitation (Yasuhara et al., 2003; Revil, 1999, 2001; Revil et al., 2006; Raj, 1982; Yasuhara and Elsworth, 2006; Elsworth and Yasuhara, 2006, 2010; Yasuhara et al., 2006b; Neretnieks, 2014). The involved mechanisms in modeling are basically that the mineral tends to dissolve where the higher stress is applied and will precipitate where the stress is lower. As a result, one may expect that such mechanisms would also contribute to the evolution of natural fractures where the two surfaces of a fracture are often in contact with each other in some asperities and form gaps elsewhere, as shown schematically in Figure 3. At contact points, the local stress is higher than the average stress on the whole fracture since stress is concentrated on a small part of the entire fracture surface. At points where the local stress is much larger than the asperity strength, mechanical crushing may occur. This will produce small mineral fragments within the fracture. When fractured rocks contain water, the minerals dissolve if the water is under-saturated. At contacting asperities with high concentrated stress the solubility of minerals is higher than at fracture voids where only the water pressure acts. The stressed minerals then tend to dissolve there and precipitate where the stress is lower. This “pressure dissolution” process at fracture contacts may lead to closure of the fracture as the stressed crystals shrink. In contrast, at fracture voids where the stress is lower than contact asperities the process of “free-face dissolution” may occur and result in local fracture widening if the water is under-saturated (Liu et al., 2006). Recent developments and investigations have shown that a relatively small change of fracture aperture would significantly influence solute transport in a rock fracture by changing the water flow rate and the diffusional exchange between the flow channel and the porous rock matrix (Koyama et al., 2008; Zhao et al., 2011). Therefore, it is worth to investigate how a fracture aperture would change under the combined effect of stress and flow. The aim is, by developing a mechanistic model that accounts for the contributions of both the matrix diffusion and the stagnant water zones, to understand and highlight the most important conditions that govern the evolution of fracture aperture in fractured rocks in the time scales relevant to the performance assessment of deep geological repositories for spent nuclear fuel. The rest of the thesis is organized as follows. In Chapters 2, the system for transport of a stable nuclide is described conceptually and formulated mathematically in the time domain. The analytical solutions in the Laplace domain are then derived for the concentration of the stable nuclide at the outlet of a flow channel. Simulations are subsequently performed to investigate the relative importance and contributions of the rock matrix layers and the stagnant water zone in retarding solute transport in fractured rocks. Chapter 3 follows the same structure as Chapter 2, but for transport of nuclides in an arbitrary length decay chain; the simulation results obtained from a rectangular channel with linear matrix diffusion are then compared with those obtained from a cylindrical channel with radial matrix diffusion. In Chapter 4, a model is presented and simulations are carried out to describe the evolution of fracture aperture under the combined effect of stress and flow. The concluding remarks are given in Chapter 5.

6

2 TRANSPORT OF A STABLE NUCLIDE

In this chapter we consider the transport of a stable nuclide in fractured rocks. The system under study is schematically shown in Figure 4, and has been well discussed in Paper I. It is an idealization of the basic building block of a network of intersecting fractures with distributions of lengths, orientations and positions (Moreno and Neretnieks, 1993b; Gylling et al., 1999). As discussed previously, the stagnant water zone may exist in the cases involving narrow channels in the fracture plane. To investigate the effect of the stagnant water zone on solute transport in fractured rocks, where the rock matrix is composed of several geological layers with different properties, a reasonable simplifying model is developed. 2Wf

2Ws

2bf

zs

z

δfi

δsi

z=0

zs = 0 x

u

2bs

y y=0

x=0

Figure 4. Flow in a channel in a fracture from which solute diffuses into rock matrix as well as into stagnant water in the fracture plane and then further into the rock matrix.

Figure 4 illustrates that the flow is assumed in the model to take place between two smooth parallel surfaces separated by a distance 2bf, forming a rectangular channel of width 2Wf. The water velocity, u, is assumed to be uniform across the flow channel. Likewise, the stagnant water zone next to the flow channel is conceptualized as a rectangular cuboid with the same length as the flow channel but different aperture and width; they are 2bs and 2Ws, respectively. The rock matrix adjacent to the flow channel is assumed to consists of nf layers with different thicknesses of 𝛿𝑓𝑖 , i = 1, 2, …, nf; the layers close to the fracture surface with i < nf constitute 7

the altered zone, while the last one represents the intact wall rock. Similarly, the rock matrix adjacent to the stagnant water zone is considered to make up of ns layers with different thicknesses of 𝛿𝑠𝑖 , i = 1, 2, …, ns; the layers close to the stagnant water zone with i < ns also constitute the altered zone, while the last one stands for the intact wall rock. When applied for practical use, however, we may simply assume that the rock matrices adjacent to both the flow channel and the stagnant water zone have identical structures and geological properties. In addition, the geological materials in the altered zone such as breccia, mylonite and cataclasite etc. may also be lumped into only one layer to share the same properties. Despite these simplifications, it should be possible to choose the flow parameters and the geological parameters such that pessimistic estimates can be obtained for the sake of safety and performance assessment of a deep repository for spent nuclear fuel. Based on the idealized geometry of the system, as shown in Figure 4, the model accounts for advection through the flow channel and diffusion of solute from the flow channel directly into the adjacent rock matrix. The solute may also at first diffuse into the stagnant water zone and then from there into the rock matrix adjacent to it. Although in Figure 4 diffusion in both the stagnant water zone and the rock matrix is shown only in the left and up directions, the model also considers diffusion in opposite directions. Moreover, linear equilibrium sorption onto both the fracture surface and the rock matrix has been included. The advection within the rock matrix is, nevertheless, ignored since the permeability of the geological layers of the rock matrix is generally low. As a consequence, solute transport in the rock matrix results solely from molecular diffusion in the model. Furthermore, a complete mixing of the solute across the flow channel has been postulated (Buckley and Loyalka, 1993) in the model, so that it is unnecessary to take the transverse dispersion in the flow channel into account. It should be noted that, when necessary, this phenomenon could be readily incorporated into the model but it was shown to have a negligible effect on solute retardation compared to the other phenomena considered. Likewise, the longitudinal dispersion is also neglected as it is generally dominated in fractured media by the differences in the residence times of solutes transported along different flow paths and hence it does not need to be treated on the level of individual flow channels (Gylling, 1997). Nevertheless, the general procedure for combining longitudinal dispersion with any retention model could be followed (Painter et al., 2008).

2.1

Mathematical model

By using the coupled 1-D approach (Hodgkinson and Maul, 1988), the system under consideration can be identified as one made up of four subsystems; the flow channel, the stagnant water zone, the rock matrix adjacent to the flow channel, and the rock matrix adjacent to the stagnant water zone. As illustrated in Figure 5, these subsystems can be characterized by nine parameter groups that describe solute transport properties within each subsystem and/or through the interface between the subsystems. The physical meanings of these characteristic parameters are given in Table 1, and a complete list of symbols, subscripts, and superscripts used in the mathematical model is given at the end of the thesis. The 1-D approach essentially models a 3-D system as an assemblage of 1-D subsystems, in which solute transport is assumed to take place only in one direction. The exchange of solute between the subsystems is, then, described as a sink or source located at the interface that

8

couples the 1-D transport equations. Since the mathematical model has been detailed in Paper I, only the governing equations are presented here for ease of reference.

Altered zone

Intact wall rock

i  Ds

MPG is

i  Df

MPG if

Fs

Ff

s

N

Stagnant water

f Flow channel

Figure 5. Schematic side view of subsystems with their characteristic parameters.

2.1.1 Solute transport in the flow channel

Based on the above considerations, the equation of continuity describing the transport of a nonradioactive solute in the flow channel can be written as, Rf

C f

C f

b D C s  u  s s t x b f W f y

 y 0

Def1 C 1pf bf

z

(1) z 0

where y and z are the coordinates directed into the stagnant water zone and into the rock matrix adjacent to the flow channel, respectively, and both are perpendicular to the channel. The subscripts f and s have been used to indicate the flow channel and the stagnant water zone, respectively. The first term at the right-hand side of Equation (1) describes the advection, while the last two terms describe the mass transfer at the interfaces between the flow channel and the stagnant water zone and between the flow channel and the first layer of the rock matrix adjacent to it, respectively. The initial condition is, C f ( x,0)  0

(2)

while the boundary condition for a time dependent concentration at the inlet is given by, C f (0, t )  Cin (t )

(3)

9

Table 1. Characteristic parameters. Notation (dimension)

Definition

Physical significance

MPGif

MPGif  Defi Ripf  ipf

The material property group, for the rock matrix layers adjacent to the flow channel, measures the ease with which solute diffuses into the pore water, but mostly the capacity to retain solute, within the matrix layer. A high value of it means that the matrix layer can hold a large amount of a given solute in both the dissolved and the adsorbed states.

MPGis  Desi Rips ips

The material property group for the rock matrix layers adjacent to the stagnant water zone has exactly the same meaning as the material property group for the rock matrix adjacent to the flow channel.

Ff (TL-1)

Ff   f RSV

The Ff-ratio is the ratio of the flow-wetted surface of the flow channel to the advection conductance (volumetric water flow rate). A high value of the Ff–ratio means that a large amount of solute carried by the flowing water can quickly be transported into the rock matrix.

Fs (TL-1)

W Fs  s Ds bs

(LT-1/2)

MPG is

(LT-1/2)

2

τf (T)

f 

τs (T)

 Dfi (T)

The characteristic time for advection in the flow channel which is actually the water residence time but can also be regarded as the solute travel time along the flow channel without accounting for any interaction with the surroundings.

x u

Ws Ds

The characteristic time for solute diffusion in the stagnant water zone which characterizes solute travel time in the direction perpendicular to the flow channel.

( if ) 2

The characteristic time for solute diffusion in the ith layer of the rock matrix adjacent to the flow channel.

2

s  i  Df 

i  Ds (T)



N (-)

N

The Fs-ratio is the ratio of the wetted surface of the rock matrix adjacent to the stagnant water zone (the stagnantwater-wetted surface) to the diffusion conductance of the stagnant water zone. It describes the competition between two important processes governing the rate of solute transport; diffusion into the stagnant water zone and diffusion from the water into the adjacent rock matrix. A large value of the Fs–ratio indicates that a large amount of solute diffusing from the flow channel into the stagnant water zone can quickly be sucked by the adjacent rock matrix.

Dafi

( si ) 2  i Das

The characteristic time for solute diffusion in the ith layer of the rock matrix adjacent to the stagnant water zone.

Ds xbs / Ws ub f W f

The N-ratio is the ratio between the characteristic rate of diffusion into the stagnant water zone and the characteristic rate of advection through the flow channel. It quantifies the fraction of solute that can depart from the flow channel into the stagnant water zone. A high value of the N–ratio indicates that the stagnant water zone has a large capability to capture the solute passing by and then deliver it from there into the adjacent rock matrix.

i Ds

10

2.1.2 Solute transport in the stagnant water zone

Assuming that the diffusion in the direction parallel to the flow is negligible, the equation of continuity describing the transport of solute in the stagnant water zone, Cs, is given by, Rs

1 C 1 C s  2 C s Des ps  Ds  2 t b  z y s s

(4) z s 0

where zs is the coordinate directed into the rock matrix adjacent to the stagnant water zone and it is perpendicular to the interface between the stagnant water zone and the adjacent rock matrix. The last term at the right-hand side of Equation (4) describes the mass transfer across the interface between the stagnant water zone and the rock matrix adjacent to it. The initial condition is, C s ( x , y ,0 )  0

(5)

and the boundary conditions are given by, Cs ( x,0, t )  C f ( x, t )

(6)

and C s y

0

(7)

y  2Ws

The processes of solute transport in the flow channel and in the stagnant water zone are thus coupled through Equation (6), which describes the continuity of the aqueous concentration of the solute.

2.1.3 Solute transport in the rock matrix adjacent to the flow channel

Accounting for molecular diffusion and linear equilibrium adsorption, the equation of continuity describing the transport of a non-radioactive solute in the pore water of the ith layer of the rock matrix adjacent to the flow channel, C ipf , is, i = 1,2,…,nf, C ipf t

 Dafi

 2 C ipf

for

z 2

where we have defined z if ,in 

z if ,in  z  z if ,out ,

i

i

k 1

k 1

(8)

  kf 1 , z if ,out    kf and  0f  0 .

The subscripts “in” and “out” indicates the beginning and the end of the matrix layer. The apparent diffusivity Dafi is defined as, i = 1,2,…,nf,

11

Dafi



Defi

(9)

R ipf  ipf

The initial conditions are, i = 1,2,…,nf, C ipf ( x, z,0)  0

z if ,in  z  z if ,out ,

for

(10)

while the boundary conditions are given, respectively, by, C ipf ( x, z, t )  C ipf1 ( x, z, t )

at

z z if ,in ,

(11)

at

z z if ,out ,

(12)

and Defi

C ipf z

 Defi 1

C ipf1 z

n 1

where we have defined C 0pf  C f ( x, t ) and Deff

0.

When i = 1, Equation (11) describes the continuity of the aqueous concentration of the solute at the interface between the flow channel and the first layer of the rock matrix adjacent to it. This couples the behavior of solute transport in the flow channel and in the adjacent rock matrix.

2.1.4 Solute transport in the rock matrix adjacent to the stagnant water zone

Similar to Equation (8), the equation of continuity describing the transport of a non-radioactive solute in the pore water of the ith layer of the rock matrix adjacent to the stagnant water zone, C ips , can be written as, i = 1,2,…,ns,

C ips t

D

i as

 2C ips z s

for

2

where we have defined z si ,in 

z si ,in  z s  z si ,out ,

i

i

k 1

k 1

(13)

  sk 1 , z si ,out    sk and  s0  0 .

The apparent diffusivity Dasi has the same definition as Dafi in Equation (9), except for changing the subscripts from f to s. The initial conditions are, i = 1,2,…,ns,

C ips ( x, z s ,0)  0

for

z si ,in  z s  z si ,out ,

(14)

while the boundary conditions are given, respectively, by,

C ips ( x, z s , t )  C ips1 ( x, z s , t )

at

z s  z si ,in ,

and 12

(15)

Desi

C ips z s

 Desi 1

C ips1

z s  z si ,out ,

at

z s

(16)

where we have defined C 0ps  Cs ( x, y, t ) and Desn 1  0 . s

When i = 1, Equation (15) describes the continuity of the aqueous concentration of the solute at the interface between the stagnant water zone and the first layer of the adjacent rock matrix. It couples, therefore, the behavior of solute transport in the two neighboring subsystems. It should be noticed that the mass transfer across the interface between the rock matrix adjacent to the flow channel and the rock matrix adjacent to the stagnant water zone is not accounted for in the model. The reason for this is that, the diffusion in the stagnant water zone is fast compared to that in the low porosity rock matrix and therefore the latter can be neglected. In other words, diffusion in the x and y directions within the rock matrix is too small compared to diffusion in the z direction.

2.2

Solution in the Laplace domain

The governing partial differential equations for solute transport given above can be solved by applying the Laplace transformation approach (Watson, 1981). It removes the time variable leaving a system of ordinary differential equations, the solutions of which yield the transform of the aqueous concentrations as a function of space variables. Since the Laplace transform technique allows one to get the aqueous concentration in the rock matrix in the Laplace domain without obtaining the concentration profiles within the flow channel and the stagnant water zone, it is preferable to begin with the Laplace-transformed equations for solute transport in the rock matrices and then continue with the equations in the stagnant water zone and the flow channel. The detailed mathematical procedure can be found in Paper I, and the results is given in what follows. In the case involving only a non-radioactive solute, the solution for the Laplace-transformed concentration in the flow channel is found to be,



 

C f ( x, s)  Cin exp   f exp  N s tanh(2 s )



(17)

where and henceforth the over-bar is used to denote the Laplace-transformed quantities, and s stands for the Laplace transform variable. In Equation (17),  f is defined as,  f  R f  f s  Ff P f s

(18)

where we have introduced three characteristic parameters τf, Ff and N. For ease of reference, these parameters together with other characteristic parameters are presented in Table 1 with their definitions and physical significance. Likewise,  s is defined for the stagnant water zone in the same way as Equation (18) except for changing the subscripts from f to s. The Pf-parameter in Equation (18) is given by, 13

Pf 

k | M kf 1 | MPG kf (1) k tanh( kf )  if cosh( if ) f | k 1 i 1 nf

 |M



(19)

where Mf is a matrix that characterizes solute transport properties within each layer of the rock matrix and through the interface between adjacent layers, and it is defined as, i 1  i i j   sinh(  ) sinh(  ) cosh( lf ) j  i  f f f l  j 1    if cosh( if ) j i  j  i MPG f j  i 1  f MPG i f  0 otherwise



 

M f  m if, j

m if, j

with

n f n f

(20) In Equation (19), Mfk+1 stands for the sub-matrix of Mf whose first to the kth rows and also columns have been deleted and |Mf| and |Mfk+1| are the determinant of Mf and Mfk+1, respectively. The first exponential function in Equation (17) describes, then, the contribution of the subsystem composed of the flow channel and the adjacent rock matrix in retarding solute transport and retention, as  f given in Equation (18) is determined solely by the parameters Rf, τf, Ff and Pf which are all related to the flow channel and the adjacent rock matrix. By the similar token, the second exponential function describes the contribution of the subsystem composed of the stagnant water zone and the adjacent rock matrix in retarding solute transport and retention with the N-ratio accounting for the fraction of solute that can depart from the flow channel into the stagnant water zone. Since Equation (17) expresses the ratio of the outlet and the inlet concentrations (in the Laplace domain) as a product of two exponential functions, it can be easily extended to the case where the transport path consists of a series of piece-wise constant channels, along each of which the surrounding media have constant geological properties. The result is, for L successive channels,  C f  Cin exp  

   if  exp  i 1   L



 N L

i 1

i



  is tanh(2  is )  

(21)

It signifies that, when applied for practical use such as in predicting the behavior of nuclide transport in fractured media from a repository to the biosphere, it is unnecessary to calculate solute concentration at every outlet of the channels by performing the inverse Laplace transform of Equation (17). Rather, one only needs to follow the flow channels along which solute travels and document the involved characteristic parameters. The information would then be sufficient to obtain solute concentration, i.e. the breakthrough curve, once at the final end of the transport path, by use of De Hoog algorithm (De Hoog et al., 1982) to numerically transform Equation (21) back to the time domain. This would not only save computation time dramatically but also give accurate predictions, provided that the assumptions underlying Equation (21) are justified to be fair for the purpose of safety assessment of a deep geological repository for spent nuclear fuel.

14

2.3

Simulations and discussion

In this section, several cases of varying complexity are simulated to investigate the relative importance and contributions of different processes and mechanisms on solute transport and retardation in fractured rocks. In the simulations, we define C f (0, t )  c0 (t ) , for a point-source boundary, where δ is the Dirac delta function, and c0 is not a concentration but the ratio between the instantaneous mass injected and the volumetric flow rate through the channel. This gives C f (0, s )  c0 , and as a result the breakthrough curve (the response function) C f c0 to the pulse injection, can readily be obtained by use of the inverse Laplace transform (Watson, 1981) of Equation (17). The data used in the base case are presented in Tables 2, 3 and 4 for the flow channel, the stagnant water zone and the adjacent rock matrix composed of five layers, respectively. The solute-dependent properties in MPG such as the diffusion coefficient and the distribution coefficient are those for cesium, while the geometry of the system and the geological properties of the rock matrix are representative of those from field observations (SKB, 2006). The database available for transport and sorption properties of the rock matrix, such as thickness and porosity, can also be found in SKB report (2003). To simplify the simulation without losing generality, however, we assume that the rock matrices adjacent to both the flow channel and the stagnant water zone have identical structures and properties.

Table 2. Flow channel properties. Symbol (unit) bf (m) Rf (–) u (m/yr) Wf (m) x (m)

Value 1.0e-04 1 5 0.1 50

Definition Rectangular channel half-aperture Retardation coefficient Groundwater velocity Rectangular channel half-width Channel length

Table 3. Stagnant water zone properties. Symbol (unit) bs (m) Ds (m2/yr) Rs (–) Ws (m)

Value 1.0e-04 0.0315 1 0.5

Definition half-aperture Diffusion coefficient Retardation coefficient half-width

Table 4. Rock matrix properties for cesium. Rock zones Intact wall rock Altered rock Cataclasite Fracture coating Fault gouge

MPG (m2/yr)0.5 0.0057 0.0138 0.0179 0.1200 0.6399

15

τD (yr) 1.43e06 3.81e04 128 0.02 0.7

A step-by-step procedure is, then, adopted to carry out sensitivity analysis, by changing one of the characteristic parameters at a time. As detailed in Paper I, it starts from a simple case where the rock matrix does not consist of any layers of altered zone with the aim of exploring the effect of the intact wall rock. The rock matrix is thereupon extended to include one or four more layers in order to facilitate understanding of the effect of the altered zone. Following this, the stagnant water zone is taken into account but assuming that the rock matrix consists only of the intact wall rock. Thereafter, the results with respect to the effects of several matrix layers and also stagnant water zone are shown. In most cases, solute transport in fractured rocks would take place in such a system that consists of not only the flow channel, the stagnant water zone, the intact wall rock, but also the altered zone which may be composed of altered rock, cataclasite, fault gouge and fracture coating etc. all at the same time. This, or even more complex, system would certainly make it hard to clarify the effects of different subsystems in retarding solute transport. However, based on sensitivity analysis made in Paper I, it is anticipated that each of the subsystems would play similar roles in determining the fate of solute transport even when everything is put together into the system. To illustrate this, we compare in Figure 6 the breakthrough curves obtained for different cases where the rock matrix consists of either none, one or four layers of altered zone, in addition to the intact wall rock, with and without accounting for the existence of the stagnant water zone. It is assumed that the altered zone consists solely of the altered rock in the case where only one matrix layer of altered zone is present; otherwise it is composed of altered rock, cataclasite, fault gouge and fracture coating.

Figure 6. Comparisons of breakthrough curves when there is/there is not stagnant water zone and rock matrix is composed of two/five layers.

It is seen in Figure 6 that when the stagnant water zone is not accounted for, the altered rock itself not only weakens the peak value of Cf (Cf,peak) but also delays the time at which the response function reaches the peak (tpeak). When other layers of altered zone are also included, however, the breakthrough curve only slightly shifts to the right, resulting in a minor decrease in Cf,peak. This indicates that the layers other than the altered rock have a moderate effect in 16

retarding solute transport, although they have large material property groups, they are too thin to have a significant influence on solute retardation. More importantly, it is noted that the effects of the stagnant water zone in determining the behavior of the breakthrough curves are always similar irrespective of how the rock matrix is structured. It always results in a considerable decrease in Cf,peak and a significant increase in tpeak, even if the altered zone consists solely of cataclasite, fault gouge or fracture coating instead of the altered rock. The reason for this is that the structure of the rock matrix only influences the equivalent material property group to some extent, if we take it as a whole as an equivalent intact wall rock, but not the N–ratio, the Fs–ratio and the diffusion time  s . It is, however, the latter characteristic parameters that play the key roles in determining the effects of the rock matrix adjacent to the stagnant water zone in retarding solute transport under the condition that the intact wall rock is thick enough to act as a deep trap for solute. Only in the cases with a wide channel and a narrow stagnant water zone so that F𝑠 MPG1𝑠 is far smaller than F𝑓 MPG𝑓1, the influence of the stagnant water zone in retarding solute transport might become marginal. The smaller the F𝑠 MPG𝑠𝑚 , the less the rock matrix adjacent to the stagnant water zone would contribute to the retardation of solute transport. Therefore, it can be concluded that the key process in retarding the solute transport is nothing but the matrix diffusion. The presence of the stagnant water zone allows the solutes to access a larger fracture surface from which diffuse into the rock matrix continues. An additional space becomes, then, available for solute transport in the rock compared to the simple system without stagnant water zone. Although the results are shown and discussed only for a stable species in this chapter, the same effects are expected to hold also for the transport of radionuclides in fractured rocks.

17

3 TRANSPORT OF NUCLIDES IN AN ARBITRARY LENGTH DECAY CHAIN

In this chapter, we consider migration of nuclides in an arbitrary length decay chain for two cases of channel geometry: rectangular flow channel with linear matrix diffusion (LMD) and cylindrical flow channel with radial matrix diffusion (RMD). The rock matrix is assumed in both cases to be composed of nf geological layers with different thicknesses of 𝛿𝑓𝑖 , i = 1, 2, …, nf, as schematically shown in Figures 7 and 8, respectively. To simplify the problem, somehow, the stagnant water zone and its adjacent rock matrix are assumed to be absent in the system studied, in order to facilitate comparison of the modeling results with those from existing conventional models. 2Wf

2bf

z

δfi

z=0 x

u

x=0

Figure 7. Flow in a flat channel in a fracture from which solute diffuses into the adjacent rock matrix layers.

Noticeably, the channel width has been found to vary typically from a few mm to several tens of cm (Moreno and Neretnieks, 1991; Abelin et al. 1983, 1985, 1994; Tsang and Neretnieks, 1998). This fact seems to favor the application of the model with RMD for practical use. However, in field experiments the observation times are usually short, which makes the impact of radial matrix diffusion on solute retardation to be small or negligible for wide channels. Therefore, the model with LMD is more applicable in these cases. The strong effect of RMD will become more pronounced mostly after tens to thousands of years, which is the time scale of interest of this work for the sake of performance assessment of deep repositories for spent 18

nuclear fuel. In this study, comparisons between the results obtained from the model with LMD and those from the model with RMD will be made, to explore the relative impact and contributions of the altered rock matrix and the decay chain on transport of radionuclides in fractured rocks.

δf3 δf2 δf1 rr0 f

Matrix Diffusion

x=0

Figure 8. Flow in a cylindrical channel from which solute diffuses into the adjacent rock matrix layers.

3.1

Mathematical model

Again, by using the coupled 1-D approach, as was done by Hodgkinson and Maul (1988) and Neretnieks (2006), two subsystems can be identified in the systems shown in Figures 7 and 8, respectively; the flow channel and the rock matrix composed of different geological layers, for which the equations of continuity describing the transport processes can be formulated individually. These two subsystems can also be characterized by those parameters given in Table 1, as detailed in Papers II and III. Without further discussions, the following sections simply present the systems of governing equations, which are valid for both models involving either rectangular or cylindrical channels.

3.1.1 Solute transport in the flow channel

Given a radioactive nuclide decay chain of length m, the equation of continuity describing the transport of the jth nuclide in the flow channel can be expressed as, j = 1,2,…,m, Rf ,j

C f , j t

 u

C f , j x

 RSV De1, j

C 1p, j r

 R f , j  j C f , j  R f , j 1 j 1C f , j 1

(22)

r rf

where we have defined λ0 = 0, and RSV as the ratio of the flow-wetted surface of the flow channel to its volume. For the rectangular channel where the flow takes place between two parallel surfaces separated by a distance 2bf (see Figure 7), the parameter RSV would become 1/bf, while

19

for the flow channel with right cylindrical geometry of radius rf (see Figure 8), the entity RSV would be equal to 2/rf. The first term at the right-hand side of Equation (22) describes the advection through the flow channel, whereas the second term indicates the mass transfer at the interface between the flow channel and the rock matrix adjacent to it. The third term represents the mass loss of member j by decay, and the forth term denotes the mass ingrowth of member j by radioactive disintegration of member j-1. Noticeably, this equation is valid for both cases involving either rectangular or cylindrical channels. In the former case, the variable r will be changed to z, and consequently the boundary r = rf will be replaced by z = 0. Equation (22) is to be solved with the initial condition that, j = 1,2,…,m,

C f , j ( x,0)  0

(23)

and the boundary condition,

C f , j (0, t )  Cin, j (t )

(24)

3.1.2 Solute transport in the rock matrix

Accounting for diffusion and linear equilibrium adsorption, the transport equation describing the transport of nuclide for the species j in the pore water of the ith layer of the rock matrix, C ip, j , is given by, j = 1,2,…,m, i = 1,2,…,nf, R ip , j

C ip , j t

D ip , j   q C ip , j r  q r r r 

   R ip , j  j C ip , j  R ip , j 1 j 1C ip , j 1  

for

i rini  r  rout

(25) i 1

i

k 0

k 0

i    k and 𝛿 0 = 𝑟𝑓 . where we have defined rini    k , rout

In Equation (25) the exponent “q” is equal to 0 in the case of rectangular channels and 1 in the 𝑖 𝑖 case of the cylindrical channels. Correspondingly, the boundaries 𝑟𝑖𝑛 and 𝑟𝑜𝑢𝑡 for the case of 𝑖 𝑖 cylindrical channels will be replaced by the boundaries 𝑧𝑓,𝑖𝑛 and 𝑧𝑓,𝑜𝑢𝑡 , respectively, in the case of rectangular channels. Equation (25) is to be solved with the initial condition that, j = 1,2,…,m, i = 1,2,…,nf,

C ip, j ( x, r,0)  0

for

i , rini  r  rout

(26)

at

r  rini

(27)

and the boundary conditions,

C ip, j ( x, r, t )  C ip,1j ( x, r, t ) and

20

Dei , j

C ip , j r

 Dei , j1

C ip,1j r

at

i r rout

(28)

n 1

where we have defined C 0p, j  C f , j ( x, t ) and De,fj  0 . Comparison of Equations (22) to (24) with Equations (1) to (3), and Equations (25) to (28) with Equations (13) to (16), shows that, when radioactive decay is accounted for, the only modification is to add one decay term and one ingrowth term into the governing equations. The initial and the boundary conditions keep invariant for the subsystems involved. This principle applies also to the case when the stagnant water zone and its adjacent rock matrix are taken into account. As mentioned before, for the moment, we only consider a simplified system that consists of a flow channel and an adjacent rock matrix. The aim is to explore the effect of different geological layers on the transport of radioactive nuclides, and to compare the difference between the radial matrix diffusion from a cylindrical flow channel and the linear matrix diffusion from a rectangular flow channel.

3.2

Solution in the Laplace domain

The governing differential equations for solute transport given in the mathematical model can be solved by applying the Laplace transformation approach (Watson, 1981), as was done also for a stable nuclide in the previous chapter. The detailed mathematical procedure can be found in Papers II and III. The final solution is presented in the following. In the case involving a radioactive decay chain, the Laplace-transformed concentration vector at the outlet of the flow channel where the effects of the stagnant water zone and its adjacent rock matrix have been neglected can be written as,



~ Cf  Θ exp T Θ 1 Cin

(29)

which can be transformed back to the time domain numerically by use of e.g. the De Hoog algorithm (De Hoog et al., 1982). The solution given in Equation (29) is valid for both rectangular and cylindrical channels. In ~ ~ this equation, the matrices Τ , Θ 1 and Θ are obtained from similarity transformations. Τ is a diagonal matrix, Θ is a lower triangular matrix, and Θ 1 , which is the inverse of Θ matrix, is also a triangular matrix. The reader is referred to Papers II and III for the definition of these matrices and the general procedure to get the Laplace-transformed concentration in the rock matrix and through the flow channel. The explicit analytical expressions for these matrices are, in general, difficult to obtain. Nevertheless, they are available in some simplified cases. For instance, in the case involving only one geological layer of the rock matrix and one single ~ nuclide, it can be shown that Θ  1 , Θ 1  1 , and the matrix Τ simplifies to, ~ T  R f,1τ f s  λ1   Ff Z11MPG11 s  1

(30)

The parameter Z11 in Equation (30) is defined, for the rectangular channel, as,

21

Z11  tanh(  1D,1 ( s  1 ) )

(31)

1 The two characteristic parameters 𝜏𝐷,1 and MPG11 follow the general definitions given in Table 𝑖 1. The characteristic diffusion time, 𝜏𝐷,𝑗 , and the material property group, MPG𝑗𝑖 , for the jth nuclide through the ith layer of the rock matrix are then defined, respectively, as,



i D, j

  

i 2 f Dai , j

(32)

and MPG ij  Dei , j R ip , j  ip , j

(33)

For the cylindrical channel, however, the parameter Z11 in Equation (30) becomes, Z11

~ ~ 1 ~ ~ 1 K1 (h11,1rin1 ) I1 (h11,1rout )  I1 (h11,1rin1 ) K1 (h11,1rout )  ~1 1 ~1 1 ~1 1 ~1 1 K 0 (h1,1rin ) I1 (h1,1rout )  I 0 (h1,1rin ) K1 (h1,1rout )

with

~ h11,1 

R1p ,1 s  1  D1p ,1 (34)

where Iα and Kα are the modified Bessel functions of the first and the second kind of order α, respectively. It follows from Equations (29) and (30) that,



C f  Cin exp  R f ,1 f s  1   Ff Z11MPG11 s  1



(35)

This equation describes the Laplace-transformed concentration of a single nuclide along a rectangular/cylindrical flow channel with the rock matrix composed of only one geological layer. The first term inside the exponential function in Equation (35) describes the contribution of advection in the flow channel to solute transport, as  f is the water residence time through the flow channel. The second term describes the contribution of diffusion and adsorption in the rock matrix to solute transport and retardation, since the entity Z11 , as given in Equations (31) and (34) for the rectangular and cylindrical channels, respectively, is related only to the rock matrix properties. It can, thus, be justified that, by performing the inverse Laplace transform of Equation (29), the outlet concentration in general cases involving not only decay chain but also different geological layers for the rock matrix would be of the following form, C f  f ( R f  f , Ff , MPG ij , Di , j ,  j , t )

(36)

22

This suggests that the characteristic parameters, as defined in Table 1, and the decay constants would be sufficient to describe the fate of nuclide transport through a discrete fracture in crystalline rocks. On the other hand, it is noted that Equation (29) can easily be extended to a more complex case, where the solutes transport through a series of channels with the properties of adjacent rock matrices varying from channel to channel. The result for the outlet concentration from L successive channels would, then, be given by, Cf 

~  Θ expTΘ  L

1

(37)

L i 1 Cin

i 1

This indicates that, in assessing nuclide transport from a repository to the biosphere, the characteristic parameters recorded for each of the channel-matrix blocks would suffice to get the solute concentration at the final end of the transport path.

3.3

Simulations and discussion

To account for the effect of radioactive decay, we have simplified the system conceptualized in Figure 4 to the one shown in Figure 7, where the stagnant water zone and its adjacent rock matrix have been ignored. As detailed on Papers II and III, the difference between the effect of the matrix diffusion of solutes from a rectangular flow channel and from a cylindrical flow channel was particularly investigated. Similarly for a stable nuclide in the previous chapter, in the simulations in this section a pulse-source boundary is considered for all nuclides at the channel inlet. Analysis of the response function, defined as C f , j c0, j , at the channel outlet would, then, give an insight into different mechanisms in retarding nuclide transport in fractured rocks.

Table 5. Rock matrix properties.

Symbol (unit) Dp (m2/s) De (m2/s) Rp (–) δ (m) εp (%)

Altered rock 9.28e-12 5.57e-14 6.66e03 0.2 0.6

Intact wall rock 6.13e-12 1.84e-14 1.33e04 1.0 0.3

Definition Pore-water diffusivity Effective diffusivity Retardation factor Thickness of layer Porosity

Table 6. Radionuclides properties. Symbol (unit) Kd (m3/kg) Rf (–) t1/2 (yr)

Am243 1.48e-02 1.0 7.37e03

Pu239 1.48e-02 1.0 2.41e04

23

Definition Volumetric sorption coefficient Surface retardation factor in fracture Half-life

To be specific, we consider a case where a two-member decay chain, Am243→Pu239, migrates through a rectangular/cylindrical flow channel with simultaneous linear/radial diffusion into the rock matrix comprising two geological layers. The geometrical data and physical properties of the flow channel and the rock matrix layers are listed in Tables 2 and 5, respectively (SKB, 2006; SKB, 1999). The half-width of the rectangular channel is Wf = 0.05 m and the radius of the cylindrical channel is rf = 2Wf/π. The nuclide properties for the two-member decay chain are given in Table 6 (SKB, 1999; SKB, 2010). The other data have been chosen to resemble the conditions in water-saturated granite rocks hosting a nuclear waste repository at about 500 m depth underground. The bulk density of the rock matrix is assumed to be 2700 kg/m3. A stepwise procedure is also used to carry out sensitivity analysis, by changing or introducing a parameter value at a time. As detailed in Papers II and III, we start from a simple case, where the rock matrix consists of only one layer (intact wall rock), and decay is not considered for the only nuclide involved, i.e. plutonium. The rock matrix is thereafter extended to have one more layer (the altered rock), to explore the effect of the geological layer on radial matrix diffusion (RMD) and also on linear matrix diffusion (LMD). Following this, the decay chain, Am243→Pu239, is taken into account to explore in detail the difference between the effect of radial and linear matrix diffusions on nuclide chain migration. For the purpose of comparison, we kept the water velocity in both channels for radial and linear matrix diffusions identical (Neretnieks, 2015) so as to have the same water residence time  f  x u in channels. In general, however, the effect of  f on retarding solute transport is negligible when the effect of matrix diffusion is large.

Matrix Diffusion

Matrix Diffusion

bf

rf

2bf

x=0 2Wf x=0

Figure 9. Rectangular channel with linear matrix diffusion, LMD (left) and cylindrical channel with radial matrix diffusion, RMD (right); r f  2W f  .

To keep the same flow-wetted surface and water velocity, the cylindrical channel can be conceptualized as concentric cylinders, as shown in Figure 9, where the inner shell (dashed line) is impermeable and the outer shell is of radius rf = 2Wf/π with Wf being the width of the corresponding rectangular channel. The entity RSV, the ratio of flow-wetted surface of the channel to its volume, in Equation (22) would then become 1 b f for both rectangular and cylindrical flow channels.

24

As discussed above, in a simple case where only one layer involved in the homogeneous rock matrix and only one single nuclide is considered, the outlet concentration is given in Equation (35), where the only difference between rectangular and cylindrical channels is in the parameter Z11 . In the case of linear matrix diffusion, it is simply a hyperbolic tangent term given in equation (31). As a result, Z11 will approach unity in the limiting case as  1D ,1   for an infinite extent of the intact wall rock. Subsequently, for a point-source boundary, the breakthrough curve (the response function), C f c0 , can readily be obtained by applying the inverse Laplace transform (Watson, 1981) on Equations (35) and (31), for t  R f ,1 f , as, Cf

 (Ff MPG 11 ) 2  c0  exp   1t   4(t  R f ,1 f )  2  (t  R f ,1 f ) 3   F f MPG 11

(38)

In the case of radial diffusion, however, the entity Z11 given in equation (34) will approach ~ ~ K1 (h11,1rin1 ) K 0 (h11,1rin1 ) for infinite thickness of the intact wall rock. This term can be approximated as, by using series expansions of the modified Bessel functions (Abramowitz and Stegun, 1968) for arguments greater than 0.2, Z11

~ K1 (h11,1rin1 ) 1  ~1 1  1  ~1 1 K 0 (h1,1rin ) 2h1,1rin

(39)

Therefore, for short time periods, i.e. when s goes to infinity, the approximate analytical solution for a cylindrical channel, can be obtained as, for t  R f ,1 f , C f c0 

 (Ff MPG 11 ) 2 De1,1 (Ff ) 2  exp    1t   4(t  R f ,1 f )  4 f 2  (t  R f ,1 f ) 3   F f MPG 11

(40)

Comparison of Equations (38) and (40) shows that the difference between linear and radial matrix diffusion for a single nuclide makes itself felt only when exp  De1,1 (Ff ) 2 4 f starts to deviate significantly from 1. This means that, when a rock matrix with large MPG is concerned, the radial matrix diffusion will retard solute transport more efficiently than the linear matrix diffusion.





On the other hand, in the limiting case as  1D,1  0 or equivalently as the penetration depth

 1D,1 (s  1 ) . As a result, the cylindrical and rectangular channels would make no difference in retarding nuclide transport and consequently we can write, by applying the inverse Laplace transform (Watson, 1981) on Equation (35), 1  1f  0 , the entity Z1 in both Equations (31) and (34) will approach

C f c0  exp( 1t ) (t  R f ,1 f  Ff MPG 11  1D,1 )

(41)

This is the response function to the pulse injection for both rectangular and cylindrical channels when the characteristic diffusion time into the rock matrix is very small, i.e.,  1D,1  0 . This analysis indicates that the difference between the two models comes mainly from the role of matrix diffusion in retarding solute transport. With small penetration depth, the role of rock 25

matrix on drawing the solutes off from the flow channel is trivial compared to the role of flowing water on carrying the solutes along the channel. As a consequence, the two models would give identical results. To better understand this, one may also conclude from Figure 9 that, when only one matrix layer (intact wall rock) is involved, the matrix volume in the linear diffusion case is 4𝑊𝑓 𝛿 1 ∆𝑥 while it becomes 4𝑊𝑓 𝛿 1 ∆𝑥 + 𝜋(𝛿 1 )2 ∆𝑥 in the radial diffusion case. Hence, when 𝑊𝑓 ≫ 𝛿 1 , or equivalently 𝜏𝐷 ⁄𝑊𝑓 → 0, the two cases would give essentially the same matrix volume accessible for nuclide transport, leading to nearly identical outlet concentrations. With large penetration depth (compared to the width Wf or rf of the channel), however, the radial matrix diffusion from a cylindrical flow channel would give solutes a much larger water volume to transport and reside than the linear matrix diffusion from a rectangular flow channel. As a result, it is not surprising to see in the simulations that much more solutes would be retained in the rock matrix when radial diffusion takes place. This is especially true when the diffusion penetration depth has become comparable to or exceed the width of the channel such that the term De1,1 (Ff ) 2 4 f in Equation (40) begins to influence the result.

Figure 10. Comparison of the response function, for Pu239, between linear matrix diffusion (LMD) and radial matrix diffusion (RMD) where the rock matrix is composed of one/two geological layers and decay chain is not/is accounted for.

To further illustrate what happens in more general cases, we shown in Figure 10 the response function for Pu23 as either a stable species or a nuclide in a first-order radioactive decay chain of Am243→Pu23. A pulse-source boundary is considered for both nuclides at the inlet. The outlet concentration of plutonium when the rock matrix is composed of either one (intact wall rock) or two layers (an altered rock zone and the intact wall rock) is shown for the both cases of LMD and RMD. It is seen, from the left figure, that the 0.2 m layer of altered rock plays an important role in both rectangular and cylindrical models not only in decreasing the outlet concentration but also in delaying the arrival time of nuclide at the channel outlet. The right figure illustrates that accounting for two matrix layers and when decay chain is considered, the outlet concentration of plutonium obtained from the cylindrical model is much lower than that obtained from the rectangular model. The difference between the output concentrations from the two cases will increase strongly by either increasing the travel distance or decreasing the water velocity in the flow channel. As 26

discussed above, increasing the travel distance or decreasing the flowing water velocity will make much more volume of the rock matrix accessible for nuclides to diffuse into and to retard there in both the dissolved and adsorbed states for a longer time. This again indicates the important role of matrix diffusion, which has a larger impact in radial diffusion since the cylindrical channel, compared to the rectangular channel, will allow a larger volume of the rock matrix to become available for retarding nuclide transport.

27

4 EVOLUTION OF FRACTURE APERTURE

In this chapter, we wish to explore how the aperture of a channel could evolve over time and along the channel by dissolution/precipitation of the minerals of the rock, and to explore how far along a channel changes in solute concentration in infiltrating water can propagate over time. For this purpose, we develop a reasonably simple model in order to understand and highlight some important conditions and phenomena governing the evolution of fracture aperture. The system under study is schematically shown in Figure 11, and has been discussed in more detail in Paper IV. A real flow path is obviously much more complex. The fracture aperture varies not only along it but also across it. Therefore, complex flow paths rather than straight flow channels are expected to form in fractures with variable apertures (Liu and Neretnieks, 2005). As a first attempt, however, we assume that the flow initially takes place between two parallel surfaces separated by 2bf, forming a rectangular channel of width 2Wf, with a constant mean velocity of u. As shown in Figure 11, the stagnant water zone next to the flow channel is similarly conceptualized as one rectangular channel of 2bs aperture and 2Ws width. Moreover, it is assumed that the channels in the fracture are held open by contact of the fracture asperities on the fracture surfaces. The minerals at the asperities are subject to rock stresses, which are considered to be present only in the stagnant water zone adjacent to the flow channel. The rock matrix is porous and hence the mineral grains there can also dissolve or act as precipitation surfaces depending on local conditions of over- or under-saturation. . In practice, it has been shown that the rock matrix is actually a strong source/sink for dissolved minerals such as quartz in the absence of flow, when the matrix porosity is larger than a few tenths of percent (Neretnieks, 2014). In addition, it was found that when the stressed crystals are small, less than a few mm, the diffusion in the thin water film between the crystal and the wall is fast, leading the concentration to be essentially equal to that surrounding the crystal. On the other hand, the volume of water film between the dissolving stressed crystal and the fracture surface is so thin that it can be neglected compared to the crystal size (Revil, 2001). The study of Neretnieks (2014) also suggested that there is no need to consider the details of the transport of dissolved minerals in this thin film because the concentration in the thin film between the crystal and the fracture plane is essentially the same everywhere and the same as the concentration in the water surrounding the crystal when there is a strong sink/source by matrix diffusion. Under these conditions, the dissolution rate of stressed faces of the crystal depends practically only on the stress and the concentration in the water just outside the thin 28

film. In the present model the stressed crystals make up a stressed surface a fraction, α. Should the crystals or the contact surface be very large, it may be necessary to account for the radial concentration gradient in the thin water film.

Cpf (t,x,-,z)

2bs

Matrix diffusion

z=0 z1 = 0 Diffusion

2bf

Lz

Cf (t,x,-,-) Ly

z

Lx

x

y

x=0 2Wf

2Ws

y=0

Dissolving free faces

Dissolving free faces Matrix diffusion

Cps(t,x,y,z1)

Dissolving free faces 2bs

Dissolving stressed faces

Cs(t,x,y,-) Cfilm(t,x,y,-) σP

side view Figure 11. Flow in a channel in a fracture where crystals are stressed in the stagnant water zone in the fracture plane. The mass dissolved by pressure dissolution, can precipitate on the unstressed surfaces and also diffuse into the adjacent porous rock matrix. Diffusion in the rock matrix and the stagnant water zone is shown only in the upward and right direction but the model considers also diffusion in opposite directions.

With these understandings, the model accounts for free-face dissolution/precipitation of minerals that make up the surfaces of the flow channel. In the stagnant water zone, pressure dissolution occurs on the stressed surfaces and free-face dissolution/precipitation on the unstressed surfaces. Dissolved minerals in flowing water and in stagnant zone can also diffuse into and out of the adjacent rock matrix. Additionally, there is mineral exchange between the flowing water and the stagnant water zone by diffusion. The effect of any hydrodynamic dispersion is, however, neglected at this stage.

29

4.1

Mathematical model

In order to devise a model that can be solved analytically, a starting point is to consider a case where the stressed area and free-face area are constant in time, or at least the rate is so slow that pseudo- steady-state conditions can be used as an approximation. This simple model will allow us to gain some insights into which processes and mechanisms will have the larger impact on evolution of fracture aperture under different circumstances. The strength of the analytical approach is that the rates of different competing processes can be summarized in a limited number of parameter groups. The analytical approach can only handle linear processes and nonlinear processes have to be handled in a different way. The analytical solution for concentration will be used in a pseudo-steady-state approach, PSSA, to study the rate of fracture closure/opening due to stress and the conditions that may lead to the growing of the flow channel aperture. We start from a case with constant aperture in the flow channel and the stagnant water zone over time. According to the coupled 1-D approach, the system shown in Figure 11 can be identified as one composed of four subsystems; the flow channel, the stagnant water zone, the rock matrix adjacent to the flow channel, and the rock matrix adjacent to the stagnant water zone. The equations of continuity for each of subsystems can be formulated individually, and be coupled at boundaries between the subsystems. The mathematical model has been detailed in Paper IV. The main governing equations are, however, presented in the following sections. The equations describing the transport of the dissolved mineral in different subsystems are similar to those describing the transport of solutes presented in Chapter 2. They differ only in additional terms indicating dissolution processes in different subsystems.

4.1.1 Transport in the flow channel

The equation describing the transport of dissolved minerals in the flow channel, C f , can be written as, C f t

 u

C f x



 1   pf

 bk C f

eq , F



Cf  

bs Dws C s b f W f y

  pf y 0

D pf C pf bf

z

(42) z 0

with    1  C eq , F  C eq ,1 expVm F  RT  

(43)

Where Ceq,1 and Ceq,F are the equilibrium concentrations of the mineral at the reference stress  1 and at the stress  F , respectively. Equation (43) shows the influence of stress on solubility (Stumm and Morgan, 1996). The stress enhancement factor is the exponential term in which Vm is the molar volume of the mineral, R and T are the gas constant and temperature, respectively.  1 is the reference stress of 1 bar (0.1 MPa) and  F is the stress applied on the free faces of the crystals which is practically the pore water pressure at the prevailing depth, i.e., 30

 F  Pp

(44)

In equation (42), x, y and z are coordinates along the flow channel, into the stagnant water zone and into the rock matrix adjacent to the flow channel, respectively.  is the porosity of the stagnant water zone. k is the mass transfer coefficient, which is the same in all four subsystems (see Appendix A in Paper IV). The term on the left side of equation (42) accounts for mineral accumulation in the water in the flow channel. On the right side, the first term accounts for transport of dissolved mineral by advection, the second term shows free-face dissolution of the fracture wall into the flow channel, the third term represents the diffusion of dissolved mineral into the stagnant water zone in the y direction and the last term describes exchange of dissolved mineral with adjacent rock matrix by diffusion in the z direction. Equation (42) is to be solved with the initial condition that C f (t  0, x)  Ceq, F

(45)

and the boundary condition that C f (t , x  0)  C f ,in

(46)

4.1.2 Transport in the stagnant water zone

Assuming that crystals in the fracture plane are rectangular cubes of size Lz and Ly = Lx and have two stressed sides parallel to the fracture surfaces, and that diffusion in the direction parallel to the flow is negligible, the equation describing the transport of dissolved minerals in the stagnant water zone, C s , is given by, C s  2 C s 1  4 1     k  Dws   1   ps k C eq , F  C s  C eq , P  C s 2   t  L b  b y y s s   1   D ps C ps   ps  bs z1 z 0











 (47)

s

with    1  Ceq, P  Ceq,1 expVm P  RT  

(48)

Where Ceq,P is the equilibrium concentration of the mineral at the stress applied on the stressed faces of the crystal,  P , given by P 

e  PP 

(49)

with the effective stress  e defined as,

31

 e   T  PP

(50)

which is basically the difference between the overburden total vertical stress  T from the weight of the rock and the pore water pressure PP (Terzaghi, 1925). In equation (49), φ is the ratio of the contact area between stressed crystals and fracture surface to the total crosssectional area normal to the applied total stress, and it varies from 0 to 1. In equation (47), the term on the left side accounts for mineral accumulation in the stagnant water zone. On the right side, the first term accounts for transport of dissolved mineral by diffusion, the second term represents the contribution of free-face dissolution/precipitation of fracture faces as well as crystal free faces into the stagnant water zone, the third term defines the contribution of pressure dissolution of stressed faces of crystals into the stagnant water zone assuming that the fracture wall is inert where the crystal is in contact, and the last term describes exchange of dissolved mineral with adjacent rock matrix by diffusion. The ratio  is considered as a fraction of fracture wetted surface in the stagnant water zone representing all lumped stressed surfaces of crystals and it is related to  by    1  W f Ws  . In the case where there is no flow channel and only stagnant water zone is present in the fracture plane, the two ratios φ and α would be identical. The initial condition is C s (t  0, y)  Ceq, F

(51)

and the boundary conditions are C s (t , y  0)  C f

(52)

and C s y

0

(53)

y  2WS

The transport processes in the flow channel and in the stagnant water zone are, thus, coupled through equation (52) describing the continuity of the aqueous concentration of dissolved minerals.

4.1.3 Transport in the rock matrix adjacent to the flow channel

The equation describing the transport of dissolved minerals in the pore water in the rock matrix adjacent to the flow channel is C pf t

 D pf

 2 C pf z 2



1   pf

 pf



ka f Ceq, F  C pf



(54)

Where Ceq,F is the equilibrium concentration of the mineral at the stress  F in the rock matrix and a f is the specific surface area of crystals in the rock matrix adjacent to the flow channel, i.e., the surface area per unit volume of crystals. 32

In equation (54), the term on the left side accounts for accumulation of dissolved mineral in the pore water in the rock matrix. On the right side, the first term represents diffusion within the rock matrix and the second term defines free-face dissolution/precipitation of crystals in the porous rock matrix adjacent to the flow channel. Pressure dissolution of crystals in the rock matrix is not included in equation (54). The reason is that the rock matrix has usually such a low porosity, on the order of a few percent, that makes the contact (stressed) area of crystals to be large and consequently a small pressure dissolution rate which will have minor effect on the change of mineral concentration in the rock matrix. Thus, for simplicity we may assume that the whole surface of crystals in the rock matrix is subject only to free-face dissolution. The initial condition is C pf (t  0, z)  Ceq, F

(55)

and the boundary conditions are C pf (t , z  0)  C f

(56)

and C pf z

0

(57)

z  f

4.1.4 Transport in the rock matrix adjacent to the stagnant water zone

Similar to equation (54), the equation describing the transport of dissolved minerals in the pore water in the rock matrix adjacent to the stagnant water is C ps t

 D ps

 2 C ps z s

2



1   ps

 ps



kas Ceq, F  C ps



(58)

The initial condition is C ps (t  0, z s )  Ceq,F

(59)

and the boundary conditions are C ps (t , z s  0)  Cs

(60)

and C ps z s

0

(61)

z s  s

33

4.2

Fracture closure rate

The set of equations of continuity along with initial and boundary conditions, as described above, can be solved by applying the Laplace transformation approach (Watson, 1981). The detailed mathematical procedure can be found in Paper IV. The analytical solution for the concentration obtained in the Laplace space is then used to study evolution of the fracture aperture in a pseudo-steady-state, PSSA, procedure. The evolution of the fracture aperture is, as discussed in Paper IV, a complicated process governed by different mechanical deformation and chemical changes in the fracture subjected to the normal stress and fluid flow. In our conceptual model, we model the crystals to be cubes with initial size Lx = Ly = Lz and have two stressed sides parallel to the fracture surfaces as illustrated in Figure 12. The two stressed surfaces have the area 2LxLy and the four unstressed surfaces have the area 4LyLz. In the lower part of Figure 12, solid lines show the datum line for aperture of both flow channel and stagnant water zone at initial time t0 and dashed lines indicate the new size of 2bf and 2bs in different regions after a certain time t1. The closure rate of fracture in the stagnant water zone will be determined by the size Lz of the shrinking (dissolving) crystal assuming that part of the fracture wall which is in contact with the crystal is inert, i.e., not dissolving (this assumption can be made when the crystal is in contact with another mineral with much lower solubility). However, we may consider a case where the fracture wall consists of the same mineral and dissolves with the rate enhanced by the stress.

2Wf

2Ws

2Ws Flow channel Stagnant water zone

Stagnant water zone

z

Lx = Ly

x 2bf (t0)

2bs (t1)

2bf (t1)

Lz

2bs (t0)

y

Figure 12. Schematic illustration of variable fracture aperture in the flow channel and the stagnant water zone due to pressure dissolution and free-face dissolution/precipitation. Full lines show the location of the fracture surface at time t0 and dashed lines at t1.

In the following examples, we assume that part of the fracture wall, which is in contact with the stressed crystal, is inert. At the same time the free-face dissolution/precipitation in the fracture voids in the stagnant water zone may open/close the stagnant water zone to some extent. The fracture closure rate in the stagnant water zone is assumed to be mainly governed by the rate of the shrinking stressed crystals. Therefore, by neglecting the very thin water film between the stressed crystals and the fracture wall, and given that PSSA is applicable (see Appendix B in Paper IV), the fracture closure rate in the stagnant water zone can be approximated by the rate of change of Lz dissolving from two sides as, 34



dbs 1 dLz   Vm k Ceq, P  C s dt 2 dt



(62)

The flow channel aperture, which determines the hydraulic conductivity, will not only close/open with the change of the aperture in the stagnant water zone at the interface, y = 0, but also simultaneously change with the free-face dissolution/precipitation of the fracture wall in the flow channel. Consequently, the aperture evolution in the flow channel can be determined by, db f dt



dbs dt



 Vm k Ceq, F  C f



(63)

y 0

There may, thus, be free-face dissolution or precipitation depending on the local conditions of saturation, and the aperture as a whole may close if the stressed crystals in the stagnant water zone dissolve faster than the free-face of crystals of the fracture wall. As a result, the aperture changes differently at different locations along the flow path. The rate of change of the aperture varies with time and can be approximately followed using the PSSA. We will, however, only consider a series of snapshots in time in this study to follow the aperture evolution.

4.3

Simulations and discussion

In this section, a series of simulations is presented and discussed to explore the effects of different processes and mechanisms on the rate of change of the fracture aperture and mineral concentration in the water in different locations. The data used in the examples are shown in Table 7. The geometrical and physical properties of the flow channel, the stagnant water zone and the rock matrix are representative of those from field observations (SKB, 2006; SKB, 1999). The pore diffusivity in the stagnant water zone modeled as a 2-D porous medium is calculated according to Archie’s law, Dws = Dw β0.6. The reacting mineral of the rock is quartz. The solubility dissolution/precipitation related data are taken from Rimstidt and Barnes (1980), and the other data have been chosen to resemble the conditions in water-saturated granite rocks, hosting a nuclear waste repository at about 500 m depth underground. The grains in the rock matrix are assumed to have the same size as the stressed mineral crystals in the stagnant water zone. It should be noted that the fracture aperture in both the flow channel (2bf) and the stagnant zone (2bs) change with time, but in a PSSA we use time intervals which are long enough for steady state concentration profile to establish. The relative rate of the aperture change is small enough compared to the rate of concentration change for the PSSA to be used. The series of simulations will give insights into the processes dominating the stress-mediated closure/opening of fractures. The results, by using the data in Table 7 and as detailed in Paper IV, show that the time to reach steady state concentrations in the rock matrix, the stagnant water zone and the flow channel is on the order of half day to a few days, and the mineral concentration reaches steady state a few cm into each of the subsystems. More importantly, it is found that the contributions from freeface dissolution and pressure dissolution to concentration build-up in both the flowing and the stagnant water zone are nearly negligible compared to the effect of the matrix diffusion, and therefore matrix diffusion is the dominant process by 99.5% compared to the other mechanisms. 35

Table 7. Data used in the simulations Symbol (unit) af = as (m2/m3) bf (m) bs (m) Ceq,1 (mol/m3) Ceq,F (mol/m3) Ceq,P (mol/m3) Dw (m2/s) Dpf = Dps (m2/s) k (m/s)

Value 6.0e04 1.0e-04 0.5e-04 0.56 0.58 1.12 1.0e-09 5.2e-10 8.54e-12

Symbol (unit) Vm (m3/mol) Wf (m) Ws (m) α (–) β (–) δf = δs (m) εpf = εpf (–) σ1 (Pa) σe (Pa)

Value 2.26e-05 0.05 0.5 0.1 0.9 0.03 0.01 0.1e06 8.1e06

Furthermore, simulations made in Paper IV suggest that the concentration of dissolved minerals in a fracture where both free-face dissolution and pressure dissolution are present will always be between Ceq,F and Ceq,P. In the absence of pressure dissolution, the concentration will finally reach Ceq,F while when stressed crystals are dissolving the concentration will go above Ceq,F but stay below Ceq,P. It might, therefore, be concluded that the free-face dissolution generally dominates over the pressure dissolution since the available area for free-face dissolution is much larger than that for pressure dissolution. This is shown schematically in Figure 13 for the concentration profile of dissolved minerals in the stagnant water zone.

Figure 13. Schematic concentration profile in the stagnant water zone.

By using the PSSA, the system under study is now used to evaluate how the fracture aperture will change in the stagnant zone and in the flow channel. The simulations account for the fact that the stress applied on stressed faces changes over time because the stressed area changes by free-face dissolution/precipitation of crystals on their unstressed surfaces. In a fracture in which aperture changes, the stress decreases considerably if either the stressed surface grows or more stressed crystals get involved. It may be noted that in the model studied in this work, all the stressed crystals in the stagnant water zone have the same size and are stressed in the same way. The stressed surface area of crystals is initially 10% of the fracture surface area in the stagnant water zone, but will change with time. The change in the stressed area will consequently increase or decrease the stress applied on the crystals.

36

The closure/opening rate of fracture can be obtained by equations given in previous sections. At the inlet of the flow channel with incoming fresh water, the free-face crystals of the fracture wall would dissolve with a rate of Vm kC f ,in which has the value of 1.12e-16 m/s (3.5e-09 m/yr). This indicates that the channel aperture would grow by 0.1 mm in about 15,000 yr by free-face dissolution, unless the fracture in the stagnant water zone closes and brings the flow channel mechanically down with it. Therefore, the flow channel aperture increases due to free-face dissolution at least near the fresh water inlet. However, further downstream where the concentration has reached equilibrium the fracture would close slowly everywhere due to the dissolution of shrinking stressed crystals in the stagnant water zone. The rate of change of the fracture aperture is shown in Figure 14. It is seen that the fracture in the stagnant water zone will completely close in about 15,000 yr, and simultaneously the fracture aperture in the flow channel will decrease from 0.2 mm to 0.1001 mm. The aperture of the flow channel will close at the rate that the stressed crystals shrink but will grow at the rate that the free-face dissolution of crystals of the fracture wall determines. The solubility of crystals in the stressed faces increases by nearly a factor of two compared to the solubility in the free-face sides as well as those crystals in the rock matrix and the fracture wall. Therefore, the fracture aperture would be mainly determined by stress-induced dissolution of crystals in the stagnant water zone. This indicates that the fracture aperture in the flow channel will essentially decrease by mechanical effects resulting from pressure dissolution of stressed crystals in the stagnant water zone.

Figure 14. Aperture evolution with Cf,in = 0 at x = 1 m and y = 1m.

The fracture aperture of the flow channel determines the hydraulic conductivity, which affects the rates of water flow and solute transport. Thus, for safety and performance assessment of radioactive waste repositories, it can be of importance to account for the change of fracture aperture, if any, in the solute transport modeling and simulations. In the examples studied, the 0.1 mm crystals in the stagnant water zone that makes up the largest part of the fracture would have dissolved in less than 20,000 yr and no further closing will take place. In the flow channel, any change in the inlet concentration would reach steady state after less than a year and after less than 0.1 m flow distance. This indicates that essentially nothing happens more than a very short distance into the flow channel in less than 20,000 yr under the conditions considered. Therefore for the assumed conditions, and for the time scales of interest for the repositories for spent nuclear fuel, the described processes have a negligible impact on the evolution of fracture 37

aperture. However, there might be conditions under which a noticeable impact on closing/opening of a fracture over a time scale of interest in crystalline rock can be found. One example could be the increase in temperature to several hundred degrees and use of much larger flow velocities such as could be of interest in geothermal heat extraction from hot dry rock. Under these or other possible conditions, the influences of crystal dissolution, matrix diffusion, and the stagnant water zone on fracture aperture need to be addressed and studied in the future.

38

5 CONCLUDING REMARKS

In this study, four models are developed. The first three models attempt to describe reactive solute transport in channels in fractured rocks, while the fourth model is dedicated to study how fracture aperture change under combined effects of stress and flow in crystalline rocks. The first model accounts for the fact that groundwater flows only in a small part of a fracture to form one or more flow channels, while at least part of the fracture is occupied by more or less stagnant water. As a result, solutes can diffuse directly from the flow channel into the adjacent heterogeneous porous rock matrix that consists of different geological layers. They may also at first diffuse into stagnant water zones and then from there into the adjacent rock matrix layers. It is, therefore, expected that both the rock matrix and the stagnant water zone will play important roles in determining the fate of solute transport in fractured rocks. The next two models take the radioactive decay chain into account, in the cases involving either linear or radial diffusion through the rock matrix composed of different geological layers. The models allow for not only an arbitrary-length decay chain but also any number and combinations of the rock matrix layers with different properties. The solutions describing the solute transport in the flow channel, which are analytical in the Laplace domain, can be numerically inverted by use of e.g. the de Hoog algorithm (de Hoog et al., 1982) to provide the solutions in the time domain. The models can, therefore, be readily implemented in both the discrete fracture network models (Dershowitz et al., 1998) and the channel network models (Gylling et al., 1999) to describe reactive transport of solutes through channels in heterogeneous fractured media consisting of an arbitrary number of rock units with piecewise constant properties. The three models have been used to simulate a series of problems of increasing complexity to understand the collaborative effect of different processes on retarding nuclide transport in fractured rocks. The results show that, if the rock matrix consists of several layers (often more porous than the intact wall rock) with different properties, an additional space for diffusion and sorption in series to the intact wall rock will be provided for radionuclides and therefore, in the case of pulse injections, the outlet concentration would significantly change not only in decreasing the peak value but also in increasing the arrival time and the peak time. This is particularly true when the altered zone is present. Furthermore, it is shown that the presence of the stagnant water zone has a major impact on retarding solute transport. The presence of the stagnant water zone and the rock matrix adjacent to it constitute an additional space for solute

39

transport in parallel to the flow channel and its adjacent rock matrix. The governing process in retarding the solutes, when accounting for stagnant water zones, is again matrix diffusion. The simulation results for the case involving a radioactive decay chain, where the effects of the stagnant water zone and its adjacent rock matrix have been neglected, are compared for two cases involving either a rectangular flow channel with LMD, or a cylindrical flow channel with RMD. The comparison suggests that the outlet concentration from the rectangular channel is practically the same as that from the cylindrical channel in the cases where the solute penetration depth is much smaller than the channel width. For narrow channels, the diffusion into the rock matrix becomes more radial by increasing the penetration depth. This leads to a larger amount of solutes to be drawn off from the flow channel, and consequently to lower the concentration at channel outlet. This indicates that RMD is more efficient in retarding solute transport than LMD when the channel width is very small compared to the solute penetration depth into the rock matrix. Thus, the two models will have the same performance when the role of matrix diffusion on retarding solute transport is negligible. The difference between the results obtained from the cylindrical channel and those from the rectangular channel increases strongly by increasing the water residence time in the flow channel. Furthermore, the simulations show that, if a thin porous altered rock layer is considered between the flow channel and the main intact rock matrix, the outlet concentration would significantly decrease, especially for the cylindrical channel case. When the rock matrix consists of an additional layer with larger porosity, an additional water volume is provided for diffusion and sorption of solutes into the rock matrix, which will in turn decrease the concentration at the outlet of the flow channel. Certainly, this additional water volume would be larger in the case involving a cylindrical channel than otherwise. In addition, it is confirmed that ignoring decay chains can cause a major overestimation of the concentration of the daughter nuclide at the flow channel outlet. This effect is again more obvious in the case of cylindrical channels, particularly when accounting for the presence of altered rock layer in the rock matrix. To summarize, the Laplace-transformed solutions obtained in the first three models require simple input parameters, and they are very efficient computationally. The results of simulated examples indicate that the solutions can easily handle the problems of reactive transport of solutes in multilayered heterogeneous systems, such as transport from a potential repository for radioactive wastes through fractured rocks. The solutions also provide a simple means for giving insight into system sensitivities to important parameters, and they are particularly useful for validating numerical codes describing reactive transport in fractured rocks. More importantly, the simulation results highlight the fact that it is necessary to account for decay chain, stagnant water zone and also rock matrix comprising at least two different geological layers (i.e. the altered rock and the intact wall rock) in safety and performance assessment of the deep repositories for spent nuclear fuel. In addition, quantitative comparisons of the results obtained from the cases of cylindrical and rectangular channels suggest that the impact of the system geometry on breakthrough curves increases noticeably as the transport distance along the flow channel and away into the rock matrix increases. The effect of the system geometry would be more pronounced for transport of a decay chain when the rock matrix contains a porous altered layer. The fourth model is developed to evaluate the evolution of fracture apertures mediated by dissolution and precipitation of a rock mineral such as quartz. The aim is to get some basic insights into how a fracture aperture changes under combined effects of stress and flow when accounting for matrix diffusion and stagnant water zones. The model includes advective flow in the fracture that can carry in or away dissolved minerals. The model also accounts for the fact that dissolved minerals in the fracture plane, in both the flow channel and the stagnant 40

water zone, can diffuse into the adjacent porous rock matrix. The analytical solution obtained in the Laplace space is then used to study evolution of the fracture aperture under combined influence of stress and flow, in a pseudo-steady-state procedure. The simulation results give insights into some important processes and mechanisms that dominate the fracture closure or opening under different circumstances. It is found that the times involved for any changes in fracture aperture are very much larger than the times needed for concentrations of dissolved minerals to reach steady state in the rock matrix, the stagnant water zone and the flow channel. This suggests that the pseudo-steady-state model can be used to assess the evolution of concentration of dissolved minerals in the rock fracture. Moreover, it is shown that diffusion into the rock matrix, which acts as a strong sink or source for dissolved minerals, clearly dominates the rate of concentration change and consequently the rate of evolution of the fracture aperture. For the condition considered, it is found that the described processes would have a significant impact on the evolution of fracture aperture, and therefore they should be considered in safety and performance assessment of repositories for spent nuclear waste. This is an important issue not only for radioactive waste repositories, but can also be used in other areas such as geological storage of CO2 or geothermal heat extraction from hot dry rock. In short, the thesis provides some insights for understanding the effects of stagnant water zones and matrix diffusion on solute transport as well as on evolution of fracture aperture in fracture rocks. However, to simplify the problems, some assumptions were made that may cause some limitations for the application of the models and open some new possibilities for further improvements. Obviously the use of the models under different conditions (by means of e.g. comparisons between the model outputs and the results of field experiments) would provide valuable information with respect to the important mechanisms and processes that govern solute transport in fractured rocks. Attempts to improve the models should, in the long run, continue as more experimental data and field observations become available. It follows that further investigations are required to get more information regarding e.g. the size distributions of stagnant water zones in fracture plane in fractured rocks. More field observations should also be made in order to better recognize the complex structure of all rock matrix layers with distinctively different geological properties. In addition, small-scale tracer tests are expected to investigate, in particular, if radial matrix diffusion has a significant impact on retarding solute transport in fractured rocks in comparison with linear matrix diffusion. The model developed for evaluation of fracture aperture should also be validated by carrying out some experiments, and field observations might be needed to give more information on saturation conditions and crystals distributions in the fractured media. Last, but not least, the overall effects of stagnant water zones, rock matrix layers and radioactive decay chain in retarding nuclide transport in the cases of both linear and radial matrix diffusion should be carefully modeled and compared with field observations and experimental results, if any, in order to aid properly for the safety and performance assessment of repositories for spent nuclear fuel.

41

NOTATIONS

Dimensions are given in terms of mass (M), length (L), time (T), moles (N) and dimensionless (–).

af

Specific surface of dissolving crystals in the rock matrix adjacent to the flow channel (L-1)

as

Specific surface of dissolving crystals in the rock matrix adjacent to the stagnant water zone (L-1)

bf

Half aperture of the flow channel (L)

bs

Half aperture of the stagnant water zone (L)

c0

Instantaneous inlet concentration into the flow channel (NL-3)

Ceq,1

Equilibrium solubility of mineral at reference stress 1 bar (NL-3)

Ceq, F

Equilibrium solubility of mineral at stress  F (NL-3)

Ceq, P

Equilibrium solubility of mineral at effective stress  P (NL-3)

Cf

Concentration in the flow channel (NL-3)

Cf,j

The jth nuclide concentration in the flow channel (NL-3)

C in

Concentration at the inlet of the flow channel (NL-3)

Cin, j

The jth nuclide concentration at the inlet of the flow channel (NL-3)

C ipf

Pore-water concentration in the ith layer of the rock matrix adjacent to the flow channel (NL-3)

42

C ips

Pore-water concentration in the ith layer of the rock matrix adjacent to the stagnant water zone (NL-3)

C ip, j

The jth nuclide pore water concentration in the ith layer of the rock matrix (NL-3)

Cs

Concentration in the stagnant water zone (NL-3)

Dafi

Apparent diffusivity in the ith layer of the rock matrix adjacent to the flow channel (L2T-1)

i Das

Apparent diffusivity in the ith layer of the rock matrix adjacent to the stagnant water zone (L2T-1)

Dai , j

Apparent diffusivity of the jth nuclide in the ith layer of the rock matrix (L2T-1)

Defi

Effective diffusivity in the ith layer of the rock matrix adjacent to the flow channel (L2T-1)

Desi

Effective diffusivity in the ith layer of the rock matrix adjacent to the stagnant water zone (L2T-1)

Dei , j

Effective diffusivity of the jth nuclide in the ith layer of the rock matrix (L2T-1)

Dip, j

Pore diffusivity of the jth nuclide in the ith layer of the rock matrix (L2T-1)

Ds

Diffusivity in the water in the stagnant water zone (L2T-1)

Dws

Pore water diffusivity in stagnant water zone (L2T-1)

Ff

Ratio of the flow-wetted surface of the flow channel to the volumetric water flow rate (TL-1)

Fs

Ratio of the stagnant-water-wetted surface to the diffusion conductance of the stagnant water zone (TL-1)

k

Mass transfer coefficient (LT-1)

L0

Initial size of cubic crystals in both z and y directions (L)

Ly

Size of cubic crystals in y direction in stagnant water zone (L)

Lz

Size of cubic crystals in z direction in stagnant water zone (L)

Mcr

Molar mass of minerals, silica (MN-1)

MPG if

Material property group of the ith layer of the rock matrix adjacent to the flow channel (LT-1/2)

43

MPG is

Material property group of the ith layer of the rock matrix adjacent to the stagnant water zone (LT-1/2)

MPGij

Material property group of the of the jth nuclide in the ith layer of the rock matrix (LT-1/2)

nf

Number of the geological layers of the rock matrix adjacent to the flow channel (–)

ns

Number of the geological layers of the rock matrix adjacent to the stagnant water zone (–)

N

Ratio between the diffusion rate into the stagnant water zone and the mass flow rate through the channel (–)

r

Distance into the rock matrix (L)

rf

Radius of the cylindrical flow channel (L)

Rf

Surface retardation coefficient in the flow channel (–)

Rf,j

Surface retardation coefficient for the jth nuclide in the flow channel (–)

Ripf

Retardation coefficient of the ith layer of the rock matrix adjacent to the flow channel (–)

Rips

Retardation coefficient of the ith layer of the rock matrix adjacent to the stagnant water zone (–)

R ip , j

Retardation coefficient for the jth nuclide in the ith layer of the rock matrix (–)

Rs

Surface retardation coefficient in the stagnant water zone (–)

RSV

The ratio of flow-wetted surface of the channel to its volume (L-1)

s

Laplace transform variable (T-1)

t

Time (T)

u

Groundwater velocity (LT-1)

Vm

Molar volume of silica (L3N-1)

Wf

Half-width of the flow channel (L)

Ws

Half-width of the stagnant water zone (L)

x

Distance along the flow direction (L)

y

Distance into the stagnant water zone (L)

z

Distance into the rock matrix adjacent to the flow channel (L)

zs

Distance into the rock matrix adjacent to the stagnant water zone (L) 44



Fraction of wetted surface in stagnant water zone demonstrating all lumped stressed surfaces (–)



Porosity of the stagnant water zone (–)

 if

Thickness of the ith layer of the rock matrix adjacent to the flow channel (L)

 si

Thickness of the ith layer of the rock matrix adjacent to the stagnant water zone (L)

 ipf

Porosity of the ith layer of the rock matrix adjacent to the flow channel (–)

 ips

Porosity of the ith layer of the rock matrix adjacent to the stagnant water zone (–)

 ip, j

Porosity of the ith layer of the rock matrix adjacent to the stagnant water zone for the jth nuclide (–)

λj

Decay coefficient for the jth nuclide (L)

ρcr

Density of crystal, silica (ML-3)

1

Reference stress (ML-1T-2)

e

Effective stress (ML-1T-2)

F

Stress applied on the crystal free faces (ML-1T-2)

P

Effective stress applied on the crystal stressed faces (ML-1T-2)

τf

Characteristic advection time (T)

 Dfi

Characteristic diffusion time through the ith layer within the rock matrix adjacent to the flow channel (T)

i  Ds

i  D, j

Characteristic diffusion time through the ith layer within the rock matrix adjacent to the stagnant water zone (T) Characteristic diffusion time for the jth nuclide in the ith layer of the rock matrix (T)

τs

Characteristic diffusion time through the stagnant water zone (T)

φ

Ratio of the contact area between stressed crystals and fracture surface to the total cross-sectional area normal to the applied total stress (–)

Subscripts cr

Refers to crystals 45

f

refers to the flow channel or the rock matrix adjacent to it

s

refers to the stagnant water zone or the rock matrix adjacent to it

p

refers to the pore water in the rock matrix

Fdiss

Refers to free-face dissolution

Pdiss

Refers to pressure dissolution1

1

Other notations can be found at the end of Papers I, II, III and IV. 46

REFERENCES

Abelin, H., J. Gidlund, L. Moreno, and Neretnieks I., 1983. Migration in a single fissure in granitic rock, 1983, Scientific basis for nuclear waste management VII, Boston, Nov 14-17, 1983. Proceedings, 239-246. Abelin, H., I. Neretnieks, S. Tunbrant, and L. Moreno, 1985. Final report on the migration in a single fracture, Experimental results and evaluation, Stripa Project Technical Report 85-103, Nuclear fuel safety project, Swedish Nuclear Fuel and Waste Managemnet Co. Stockholm. Abelin H., L. Birgersson, L. Moreno, H. Widén, T. Ågren, and I. Neretnieks, 1991. A Large Scale Flow and Tracer Experiment in Granite II. Results and interpretation. Water Resources Research, 27(12), 3119-3135. Abelin H., L. Birgersson, H. Widén, T. Ågren, L. Moreno, and I. Neretnieks, 1994. Channeling experiments in crystalline fractured rocks. J. Contam Hydrol 15(3), 129–158. Barten, W., 1996. Linear response concept combining advection and limited rock matrix in a fracture network transport model, Water Resour. Res., 32(11), 3285– 3296. Beeler, N. M., and S. H. Hickman (2004. Stress-induced, time-dependent fracture closure at hydrothermal conditions. J. Geophys. Res., Vol. 109, B02211. Birgersson L., H. Widén, T. Ågren, and I. Neretnieks, 1992. Stripa Project: tracer migration experiments in the Stripa mine 1980–1991. SKB/OECD-NEA Technical Report 92-25, SKB, Stockholm. Birgersson, L., L. Moreno, I. Neretnieks, H. Widén, and T. Ågren, 1993. A tracer migration experiment in a small fracture zone in granite, Water Resour. Res., 29(12), 3867-3878. Buckley, R. L., S. K. Loyalty, 1993. Radionuclide migration in fractured media: numerical studies of advection/diffusion in a fracture and diffusion in the surrounding rock matrix. Ann. Nucl. Energy 20 (10), 701-718. Carrera. J., X.Sánchez-Vila, I. Benet, A. Medina, G.Galarza, J. and Guimerà, 1998. On matrix diffusion: formulations, solution methods and qualitative effects, Hydrogeology Journal (1998) 6:178–190.

47

Chopra, M., R. N. Nair, F, Sunny, and D. N. Sharma, 2015. Migration of radionuclides from a high-level radioactive waste repository in deep geological formations, Environ Earth Sci, 73:1757–1768, doi: 10.1007/s12665-014-3525-x. Cihan, A., and J. S. Tyner, 2011. 2-D radial analytical solutions for solute transport in a dualporosity medium, Water Resour. Res., 47, W04507, doi:10.1029/2009WR008969. Cvetkovic, V., 2010. Significance of fracture rim zone heterogeneity for tracer transport in crystalline rock, Water Resour. Res., 46, W03504, doi:10.1029/2009WR007755. De Hoog, F. R., J. H. Knight, and A. N. Stokes, 1982. An improved method for numerical inversion of Laplace transforms. SIAM J. Sci. Stat. Comput. 3 (3), 357–366. Davis, G. B., C. D. Johnston, 1984. Comment on “Contaminant transport in fractured porous media: analytical solution for a system of parallel fractures” by E. A. Sudicky and E. O. Frind. Water Resources Res. 20 (9), 1321-1322. Deng, H., Z. Dai, A. Wolfsberg, Z. Lu, M. Ye, and P. Reimus, 2010. Upscaling of reactive mass transport in fractured rocks with multimodal reactive mineral facies, Water Resour. Res., 46, W06501, doi:10.1029/2009WR008363. Dershowitz, W. S., G. Lee, J. Geier, T. Foxford, P. LaPointe, and A. Thomas, 1998. FracMan version 2.6: Interactive discrete feature data analysis, geometric, modeling, and exploration simulation—User documentation, Rep. 923-1089, Golder Assoc., Inc., Seattle, Wash. Dershowitz, W., Follin, S., Eiben, T., Andersson, J., 1999. SR 97Ð alternative models project: discrete fracture network modelling for performance assessment of Aberg, SKB Report R-9943. Swedish Nuclear Fuel and Waste Management Co., Stockholm, Sweden. Elsworth, D., and H. Yasuhara, 2006. Short-timescale chemo-mechanical effects and their influence on the transport properties of fractured rock, Pure Appl. Geophys. 163: 2051–2070. Elsworth, D., and H. Yasuhara, 2010. Mechanical and transport constitutive models for fractures subject to dissolution and precipitation, Int. J. Numer. Anal. Meth. Geomech. 34: 533– 549. Fourest, B., J. Perrone, P. Tarapcik, and E. Giffaut, 2004. The hydrolysis of protactinium (V) studied by capillary diffusion. Journal of solution chemistry 33 (8). Golub, G. H. and Van Loan, C. F., 1996. Matrix Computations, 3rd ed. Baltimore, MD: Johns Hopkins University Press, pp. 311. Grisak, G. E., and J. F. Pickens, 1980. Solute transport through fractured porous media, 1, The effect of matrix diffusion. Water Resources Res. 16 (4), 719-730. Gylling, B., 1997. Development and applications of the channel network model for simulations of flow and solute transport in fractured rock, Ph.D. thesis, Dept. Chemical Engineering and Technology, Royal Inst. of Technology, Stockholm. Gylling, B., L. Moreno, and I. Neretnieks, 1998. Simulation of Radionuclide Release from a Repository to the Biosphere: Using a Model-Coupling Concept, Nucl. Technol., 122(1), 93103. 48

Gylling, B., L. Moreno, and I. Neretnieks, 1999. The Channel Network Model - A tool for transport simulation in fractured media. Groundwater 37, 367-375. Hartley, L. J., 1998. NAPSAC technical summary document, AEA-D&R-0271, AEA Technol., Harwell, U. K. Hartley L., and D. Roberts, 2012. Summary of discrete fracture network modeling as applied to hydrogeology of the Forsmark and Laxemar sites, SKB report R-12-04 (www.skb.se). Hodgkinson, D. P., and P. R. Maul, 1988. 1-D modelling of radionuclide migration through permeable and fractured rock for arbitrary length decay chain using numerical inversion of Laplace transforms. Ann. Nucl. Energy 15 (4), 175–189. Idemitsu, K., H. Furuya, K. Murayama, and Y. Inagaki, 1992. Diffusivity of uranium in water saturated inada granite. Mat.Res. Symp. Proc. 257. Ivanov, V. V. and I. B. Popov, 2009. Diffusion of americium in aqueous solutions. Original Russian text. Jing, L., K. B. Min, A. Baghbanan, and Z. Zhao, 2013. Understanding coupled stress, flow and transport processes in fractured rocks, Geosystem Engineering, 16:1, 2-25, DOI: 0.1080/12269328.2013.780709. Koyama, T., Li, B., Jiang Y., Jing L., 2008. Numerical simulations for the effects of normal loading on particle transport in rock fractures during shear. Int J Rock Mech Min Sci; 45:1403– 9. Kumagai, M., M. Sazarashi, K. Idemitsu, A. Kagawa, H. Umeki and K. Ishigure, 2009. Diffusion of plutonium in compacted bentonites. Li SH, Chiou SL, 1993. Radionuclide migration in fractured porous rock: analytical solution for a kinetic solubility limited dissolution model. Nucl Technol 104:258–271. Liu, L., and I. Neretnieks, 2005. Analysis of fluid flow and solute transport in a fracture intersection a canister with variable aperture fractures and arbitrary intersection angles, Nucl. Technol., 150(2), 132–144. Liu, J., J. Sheng, A. Polak, D. Elsworth, H. Yasuhara, and A. Grader, 2006. A fully-coupled hydrological–mechanical–chemical model for fracture sealing and preferential opening, Int. J. Rock Mech. Min. Sci., 43: 23–36. Löfgren, M., and M. Sidborn, 2010a. Statistical analysis of results from the quantitative mapping of fractureminerals in Forsmark— Site descriptive modelling – complementary studies, SKB Rep. R-09-30. Swedish Nuclear Fuel and Waste Management Company, Stockholm. Löfgren, M., and M. Sidborn, 2010b. Statistical analysis of results from the quantitative mapping of fracture minerals in Laxemar. Site descriptive modelling—Complementary studies. SKB Rep. R-09-31, Swedish Nuclear Fuel and Waste Management Company, Stockholm. Molz, F., 2015. Advection, dispersion and confusion. Groundwater 53, doi:10.111/gwat.12338. Moreno L., and I. Neretnieks, 1991. Fluid and solute transport in a network of channels, SKB Technical report TR-91-44 (www.skb.se). 49

Moreno L., and I. Neretnieks, 1993a. Fluid flow and solute transport in a network of channels. J. Contaminant Hydrology 14, 163- 192. Moreno, L., and I. Neretnieks, 1993b. Flow and nuclide transport in fractured media: The importance of the flow wetted surface for radionuclide migration, in Chemistry and Migration of Actinides and Fission Products, edited by J. I. Kim and G. de Marsily, J. Contam. Hydrol., 13, 49– 71. Moreno, L., Gylling, B., Neretnieks, I., 1997. Solute transport in fractured media-the important mechanisms for performance assessment. J. Contam. Hydrol. 25, 283–298 Moreno, L., and J. Crawford, 2009. Can we use tracer tests to obtain data for performance assessment of repositories for nuclear waste?, Hydrogeol. J., 17, 1067–1080. Neretnieks, I., 1980. Diffusion in the rock matrix: An important factor in radionuclide retardation? , J. Geophys R. 85 (B8), 4379-4397. Neretnieks I., 1983. A note on fracture flow dispersion mechanisms in the ground. Water Resources Res. 19, 2, 364-370. Neretnieks, I., 1993. Solute transport in fractured rock-application to radionuclide waste repositories, Flow and contaminant transport in fractured rock. Academic Press, Inc. Neretnieks, I., 2006. Channeling with diffusion into stagnant water and into a matrix in series, Water Resources. Res. 42. Neretnieks, I., 2014. Stress-mediated closing of fractures—Impact of matrix diffusion, J. Geophys. Res. Solid Earth, 119, doi:10.1002/2013JB010645. Neretnieks, I., 2015. Solute transport in channel networks with radial diffusion from channels in a porous rock matrix, SKB Rep. R-15-02. Swedish Nuclear Fuel and Waste Management Company, Stockholm (www.skb.se). Painter, S., V. Cvetkovic, J. Mancillas, and O. Pensado, 2008. Time domain particle tracking methods for simulating transport with retention and first-order transformation, Water Resour. Res., 44, W01406, doi:10.1029/2007WR005944. Park, J. B., Y. Hwang, and K. J. Lee, 2001. Analytic solutions of radionuclide transport with the limited diffusion from the fracture into a porous rock matrix, Ann. Nucl. Energy, 28, 9931011. Piqué, A., F. Grandia, C. Sena, D. Arcos, J. Molinero, L. Duro, and J. Bruno, 2010. Conceptual and numerical modelling of radionuclide transport in near-surface systems at Forsmark. SRSite Biosphere, SKB Rep. R-10-30. Swedish Nuclear Fuel and Waste Management Company, Stockholm. Polak, A., D. Elsworth, H. Yasuhara, A. S. Grader, and P. M. Halleck, 2003. Permeability reduction of a natural fracture under net dissolution by hydrothermal fluids. Geophys Res Lett. 30:2020. Poteri, A., and M. Laitinen, 1999. Site-to-canister scale flow and transport at the Hästholmen, Kivetty, Olkiluoto and Romuvaara sites, Rep. 99– 15, POSIVA, Olkiluoto, Eurajoki, Finland.

50

Raj, R., 1982. Creep in polycrystalline aggregates by matter transport through a liquid phase, J. Geophys. Res., 87(B6): 4731–4739, doi:10.1029/JB087iB06p04731. Rasmuson A., and I. Neretnieks, 1986. Radionuclide transport in fast channels in crystalline rock. Water Resources Res. 22, p 1247-1256 Rahman, M. M., R. Liedl, and P. Grathwohl, 2004. Sorption kinetics during macropore transport of organic contaminants in soils: Laboratory experiments and analytical modeling. Water Resour. Res., 40, W01503, doi:10.1029/2002WR001946. Revil, A., 1999. Pervasive pressure-solution transfer: A poro-visco-plastic model, Geophys. Res. Lett. 26(2): 255–258, doi:10.1029/1998GL900268. Revil, A., 2001. Pervasive pressure solution transfer in a quartz sand, J. Geophys. Res., 106(B5): 8665–8686, doi:10.1029/2000JB900465. Revil, A., P. Leroy, A. Ghorbani, N. Florsch, and A. R. Niemeijer, 2006. Compaction of quartz sands by pressure solution using a Cole-Cole distribution of relaxation times, J. Geophys. Res., 111, B09205, doi:10.1029/2005JB004151. Rimstidt, J. D., and H. L. Barnes, 1980. The kinetics of silica-water reactions, Geochim, Cosmochim. Acta, 44: 1683–1699. Selnert, E., J. Byegård, H. Widestrand, S. Carlsten, C. Döse, and E.-L. Tullborg, 2009. Bedrock transport properties. Data evaluation and retardation model. Site descriptive modelling SDMSite Laxemar, SKB Rep. R-08-100. Swedish Nuclear Fuel and Waste Management Company, Stockholm. Selroos J-O., D.D. Walker, A. Ström, B. Gylling, and S. Follin, 2002. Comparison of alternative modeling approaches for groundwater flow in fractured rock, Journal of hydrology, 257, 174188. Sharifi Haddad, A., Hassanzadeh, H., Abedi, J., Chen, Z., 2013. Lumped mass transfer coefficient for divergent radial solute transport in fractured aquifers, J Hydrol. v 495. Stober I., and K. Bucher, 2006. Hydraulic properties of the crystalline basement, Hydrogeology Journal 15, 213–224. Stumm, W., and J. J. Morgan, 1996. Aquatic Chemistry, John Wiley. New York. Sudicky, E. A., E. O. Frind, 1982. Contaminant transport in fractured porous media: Analytical solutions for a system of parallel fractures. Water Resources. Res. 18 (6), 1634-1642. Sudicky, E. A., E. O. Frind, 1984. Contaminant transport in fractured porous media: Analytical solution for a two-member decay chain in a single fracture. Water Resources. Res. 20 (7), 10211029. Swedish Nuclear Fuel and Waste Management Company (SKB), 1999. Radio nuclide transport calculations, SKB Tech. Rep. TR-99-23. Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. Swedish Nuclear Fuel and Waste Management Company (SKB), 2002a. Final report of the TRUE Block Scale project, 1. Characterisation and model development, SKB Tech. Rep. TR02-13, Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. 51

Swedish Nuclear Fuel and Waste Management Company (SKB), 2003. Äspö hard rock laboratory, SKB Tech. Int. Progr. Rep. IRP-03-13. Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. Swedish Nuclear Fuel and Waste Management Company (SKB), 2003. Final report of the TRUE Block Scale project, 4. Synthesis of flow, transport and retention in the block scale, SKB Tech. Rep. TR-02-16, Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. Swedish Nuclear Fuel and Waste Management Company (SKB), 2006. Äspö hard rock laboratory, Äspö task force on modelling of groundwater flow and transport of solutes, Modelling of task 6D, 6E and 6F, using CHAN3D, SKB Int. Progr. Rep. IPR-06-19, Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. Swedish Nuclear Fuel and Waste Management Company (SKB), 2008. Bedrock transport properties Data evaluation and retardation model, SKB Rep. R-08-98. Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. Swedish Nuclear Fuel and Waste Management Company (SKB), 2010. Bedrock Kd data and uncertainty assessment for application in SR-Site geosphere transport calculations, SKB Rep. R-10-48. Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. Swedish Nuclear Fuel and Waste Management Company (SKB), 2011. Long-term safety for the final repository for spent nuclear fuel at Forsmark, SKB Tech. Rep. TR-11-01, Swedish Nuclear Fuel and Waste Management Company, Stockholm, Sweden. Tang, D. H., E. O. Frind, and E. A. Sudicky, 1981. Contaminant transport in fractured porous media: Analytical solution for a single fracture. Water Resources. Res. 7 (3), 555-565. Terzaghi, C., 1925. Principles of Soil Mechanics. Engineering News-Record, 95(19-27). Tsang, C.-F., I. Neretnieks, and Y. Tsang, 2015. Hydrologic issues associated with nuclear waste repositories, Water Resour. Res., 51, doi:10.1002/ 2015WR017641. Tsang, C.-F., and C. Doughty, 2003. A particle-tracking approach to simulating transport in a complex fracture, Water Resour. Res., 39(7), 1174, doi:10.1029/2002WR001614. Tsang C.-F., and I. Neretnieks, 1998. Flow channeling in heterogeneous fractured rocks, Reviews of Geophysics, 26(2), 275-298. van Genuchten, M. T., D. H. Tang, and R. Guennelon, 1984. Some exact solutions for solute transport through soils containing large cylindrical macropores. Water Resour. Res., 20(3), 335–346, doi:10.1029/WR020i003p00335. van Genuchten, M.Th., 1985, A general approach for modelling solute transport in structured soils, Proc. 17th international congress, hydrology of rocks of low permeability. Watson, E. J., 1981. Laplace Transform and Applications, Van Nostrand Reinhold Co., London. Yasuhara, H., D. Elsworth, and A. Polak, 2003. A mechanistic model for compaction of granular aggregates moderated by pressure solution, J. Geophys. Res. 108(B11), 2530, doi:10.1029/2003JB002536. 52

Yasuhara, H., and D. Elsworth, 2006. A numerical model simulating reactive transport and evolution of fracture permeability, Int. J. Numer. Anal. Meth. Geomech., 30: 1039–1062. Yasuhara, H., A. Polak, Y. Mitani, A. Grader, P. Halleck, and D. Elsworth, 2006a. Evolution of fracture permeability through fluid-rock reaction under hydrothermal conditions. Earth Planet Sci. Lett. 244:186–200. Yasuhara, H., D. Elsworth, A. Polak, J. Liu, A. Grader, and P. Halleck, 2006b. Spontaneous switching between permeability enhancement and degradation in fractures in carbonate: lumped parameter representation of mechanically- and chemically- medicated dissolution. Transp Porous Media. 65:385–409. Yu, Y. W., C. H. Chung, and C. L. Kim, 1994. Analytical solutions for a three-member decay chain of radionuclides transport in a single fractured porous rock. Journal of the Korean Nuclear Society 26 (4). Zhao, Z., L. Jing, I. Neretnieks, L. Moreno, 2011. Analytical solution of coupled stress-flowtransport processes in a single fracture. Computers and Geosciences. 37: 1437-1449. Zhao, Z., L. Liu, I. Neretnieks, and L. Jing, 2014. Solute transport in a single fracture: Impacted by chemically mediated changes, Int. J. Rock Mech. Mining Sci., 66: 69–75. Zimmerman, W., Bodvarsson, G.S., 1995. Effective block size for imbibition or absorbtion in dual-porosity media. Geophysical Research Letter, 22(11).

53