Radionuclide transport in fractured rock

Progress in Nuclear Energy. Vol. 30. No. 3, pp. 243-253. 1996 CopyrightQ 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0149-1970/96 $32.00

Pergamon 014~1970(95)00086-0

RADIONUCLIDE

TRANSPORT

IN FRACTURED

ROCK

M. M. R. Williams Electrowatt Engineering Services (U.K.) Ltd. North Street, Horsham, West Sussex, U.K.

ABSTRACT The mathematical models currently employed to study the migration of radionuclides in fractured rock are explained and compared. Some of the advantages and disadvantages of these are highlighted and a new approach based on an analogy with neutron transport is introduced.

I. INTRODUCTION Nuclear power reactors produce useful energy in the form of heat and electricity. They also produce hazardous waste in the form of fission products which may remain dangerous to life for many thousands of years. The problems associated with the disposal of this waste is an international one and a variety of techniques have been proposed for solving it. The basic difficulty is that of devising a containment system that will isolate the nuclear waste from the biosphere for a period of time that is long compared with the associated half-lives. In the past ten to fifteen years, the two most serious proposals for disposal have been burial at sea and burial in deep rock formations. Burial at sea has met with strong environmental opposition and the idea has essentially been abandoned. While deep rock disposal is not entirely free of antagonists, it seems to be the only feasible solution and most of the current effort is directed towards understanding the technical problems and costs which will arise. The main criterion in the selection of a site for a repository is that the host rock must have a history of stability, i.e. the rock formation shall have been quiescent in terms of earthquake or violent fracturing for hundreds of thousands of years. While the repository itself may occupy only some 3 or 4 km2, the surrounding rock formation may well have an area1 extent of 300,000 km2. The area should also be free of any ore-bearing rocks likely to be of interest to man, since future generations may mine these areas when all knowledge of the repository has been lost. In addition to the above requirements, there is another one of considerable importance. Namely, that the influence of groundwater be minimised. Rock formations, by their very nature, contain fissures and fractures due to the effect of forces in the earth’s crust. Inevitably, over a period of time, water will find its way into these void spaces and lead to a slow but finite advective motion around the repository. Now, while the canisters in which the waste is stored may be strong and backfilled with material of low permeability, over several thousand years the canister contents will be leached out into the local groundwater. It is therefore important to be able to predict the motion of the groundwater and its likely variation with time over perhaps one million years. 243

244

M. M. R. Wilhms

It will be noted from this preamble that the factors to be included in an assessment of a repository are many and varied. Geology of the site and human behaviour are but two of a host of sub-problems. In this paper we will concentrate on one particular aspect which has relevance to transport theory in a generalised sense, namely the modelling techniques that are employed to study the movement of radionuclides through fractured rock. Such a material may be regarded as a spatially random medium and so the problem becomes one of modelling groundwater, and the dissolved radionuclides, in a spatially stochastic system. We will describe some of the methods currently employed to understand the problems and introduce some new ideas that have proved useful. Where possible, we will try to relate the methods used to those of conventional reactor physics, thereby assisting the newcomer to the field to appreciate the inter-disciplinary nature of this area of activity. In order to put the problem into perspective, consider Fig 1 which is a possible setting for a nuclear waste repository.

Fig 1. Schematic diagram of a typical waste repository

No matter where the waste is stored, there will always be water present. Even in very deep repositories, groundwater will seep in, albeit very slowly. However, since we are talking about modelling events over a period of several hundred thousand years, even a low flow rate can lead to significant transport. Imagine, therefore, radioactive material dissolving into the surrounding water and then being carried in solution through the rock matrix. This will consist of a porous medium with cracks and faults of varying length and width as illustrated in Fig 2. We have in fact a spatially random medium, with possibly several scale lengths, corresponding to a few centimetres for porous dispersion, up to several kilometres for the fractures. The problem therefore is to model the transport of the water borne radionuclides through this medium. In fact, the problem is rather complex because, due to adsorption, some of the radioactive material will adhere to the rock matrix surface and come into quasiequilibrium with that in solution. This has the beneficial effect of reducing the effective velocity of the solute but it makes the mathematics more difficult. In the next section we shall explore the various methods used to model the behaviour of solute as it moves through the rock.

Radionuclide

transport

245

Microftrsures wth diffusion only

Fig 2. Two-dimensional view of microfissures in granite showing typical sizes of grains and fissures II MODELS OF TRANSPORT IN FRACTURED ROCK There are many different approaches to modelling radionuclide transport, but essentially each can be classified as one of three standard procedures; (1) classical advective-dispersion theory, (2) double porosity methods and (3) Monte Carlo generation of random fractures and associated tracking of the solute. More recently, a method based upon an analogy with neutron transport has been developed which shows some promise (Williams, 1992). We briefly describe the principles behind these methods below. A. Classical advective-dispersion theory The standard approach for dealing with this problem is to assume that the radionuclide concentration obeys a diffusion equation with advection, viz:

+mc-V.(“C)-AC

(1)

where C is the concentration of radionuclide at position r and time t. D is the dispersion coefficient, U the local Darcy velocity and h the radioactive decay constant. For chains of nuclides, this equation becomes a coupled set but nothing new is added and for simplicity we shall consider only the single nuclide case. Eqn(1) hides a number of difficulties. Recall that the medium through which the waterborne contaminant is moving is porous and contains various faults and fractures. Thus the socalled dispersion coefficient is actually an average measure of the ease of transport by random motion through the interstices of the rock. In fact, D is tensorial in nature, because the medium is stratified. Thus in practice, D=(DL,DT,@) where DL is a longitudinal component and DT is a transverse one. This is not unlike the diffusion asymmetry in a reactor due to the presence of cylindrical fuel rods. D is also found to depend, to a first approximation, linearly on the Darcy velocity. Thus we can write D=aU, where a is the so-called dispersivity or dispersion length. For example in basalt in the Hanford region, aL= 30 m and aT = 20 m.

246

M. M. R. Williams

As well as the equation for the concentration, additional relations for the Darcy velocity are needed. Now according to Darcy’s law, we have U =-K.Vh

(2)

where K is the hydraulic conductivity tensor and h is the hydraulic head z+p/pg. z being the height of water above the point of interest. In cases where V. U = 0, we have the equation V.K.Vh=O

(3)

which defines a boundary value problem for the head. This must be solved first to find U before the radionuclide equation for C can be solved. In practice, U is often prescribed based on local experimental knowledge. In this paper, U will always be assumed to be known. Eqn(1) must be modified in practical applications, and we now consider how this is done and why. Consider Fig 3, and note that waste leaches out of the container and moves along a typical fracture, crack or pore. We assume that the contaminant, in this case the radionuclides, are distributed in the water, and in the rock. The concentration of the solute in water is C (mass per unit vol of water) and in the rock it is F (mass per unit volume of rock). There are then two balance equations; one for the rock and one for the water. For the water: e$=W.[DVC-UC]-dC-I(F,C)

(4)

where E is the porosity which is constant in a saturated medium. For the rock:

(5)

(I-E)$-(I--E)AF+f(F,C)

Figure 3 Illustration of the sorption process

Radionuclide transport

241

In the rock there is no mass motion or diffusion of the contaminant. Each term is fairly self-evident, except perhaps f. This term denotes the rate at which solute is adsorbed by the rock matrix. A number of so-called ‘isotherms’ have been proposed for this function but the simplest is the linear isotherm which takes the form: f =+@C-F)

(6)

where z is a time constant. Thus there are two coupled equations to solve not unlike reactor kinetics with one group of delayed neutrons. The quantity of interest is the total concentration of radionuclide at a given position and time in liquid and rock, i.e. EC+ (I- E)F. Thus, if we add the two equations, we find: -++(I--E)F]=&V.[DVC-UC]-A[&C+(l-&)F]

(7)

In practice, the time constant z is generally very short compared with times of interest in the system, i.e. minutes, hours and days compared with hundreds of years. So it is assumed that equilibrium exists between the liquid and solid phases, i.e. F = PC. In that case the equation becomes dC

&

at=&+(l-&)/3

V.[DVC - UC] - AC

The term 1 E+(lTE)jj

=1

(1-E) +yP

1 =R

(9)

where R is called the ‘retardation factor’ and so we have finally dC at

-=

VD .

[

RVC-$k1 -AC

(10)

Thus the effective diffusion coefficient and effective advective velocity is reduced by a factor R compared with the actual ones. For some radionuclides R - 104 and so transport times are reduced considerably. Other radionuclides, however, like iodine, undergo very little retardation and R is close to unity. There is therefore a separation of isotopes during their migration through the rock. This is of course the phenomenon that accounts for the separation effect to be seen in chromatography. Thus the equation to be studied is the conventional advective - diffusion equation but with effective parameters. It should also be noted that Brenner (1991) has shown that in a chemically reactive solute, i.e. one which reacts with the host rock, the solute velocity may be greater than that of the fluid, thereby leading to a value of R which is less than unity. Such matters have received little attention. The use of such a dispersion-advection equation would seem to be not unreasonable and it forms the basis of many proprietary computer codes used to predict radionuclide behaviour over long periods of time. In order to obtain the data for this equation, experiments are performed on samples of the material and values for aT and aL deduced by least squares fitting. During the past 15 years, however, field experiments have shown that the longitudinal dispersivity aL is “scale dependent”. That is to say, one gets different values according to the

248

M. M. R. Williams

size of the experiment. Roughly speaking, the dispersivity seems to be proportional to the distance travelled. This is a most unwelcome result since it means that what was thought to be a unique property of the system, now depends on the problem. The rule of thumb for people who use the codes is to set the dispersivity equal to the length of the system under investigation! Deep faith is needed for this. There have been a number of attempts to resolve this problem and most are based upon re-interpreting the equation for the concentration as a stochastic diferentiul equation with the advective velocity U being a stochastic parameter. We have met this type of problem before in other fields. In particular, in the transmission of waves through media with random refractive index, and also in the problem of neutron transport through materials with spatially random density. To illustrate the problem, consider a two dimensional situation such that:

(11) This equation describes a medium in which there is a net flow in the x-direction with velocity U and dispersion in the x and z directions. Randomness is caused by random stratification in the z-direction which leads to a random hydraulic conductivity and hence flow velocity. Using perturbation theory, such that: C(x,z,t)= ~(x,t)+c(x,z,t) U(z) = g + u(z) D(z) = b + d(z) where c, u and d are small and random, with zero mean, it is possible to re-write the equation as:

ac

5‘

d2c -_

- d2C

a2F

at- acy 4~=4-g-u~

a??

(12)

where 5 = x-Vt is a reduced variable. This is a relatively simple stochastic differential equation with the randomness arising from an inhomogeneous source term on the right hand side. Given the statistics of the hydraulic conductivity (which affects dL and u) it is relatively easy to calculate the mean and variance in the concentration. In fact it may be shown that the equation for the average concentration can be written to a fust approximation as: (13) where A(t) is related to the fluctuations in the hydraulic conductivity. Thus the classical dispersitivity aL is augmented by a parameter A(t) which depends on time. We have a nonFickian situation and to some extent this explains the experimental observations of a varying dispersion length but is not particularly helpful in practice for quantitative assessments. B Double porosity method The classical method in section A is not unreasonable if the material is dominantly porous with little anisotropy in the dispersion. However, in practice, fractures will give rise to short cuts through the rock, i.e. fast legs, in which contaminant can can rapidly spread. The classical method in the form described above, is not suitable for such problems. To overcome this, some authors have divided the rock into a number of ideal fractures which take the form of infinite slits of finite width (Sudicky and Frind, 1982). Fig 4 shows an example.

Radionuclide

“F

uanspon

C, fracture

Fig 4. Double

249

)_z

porosity model

Equations are then constructed for the fracture and porous matrix separately. The equation for the fracture is written

acF _ DF a2cF at-R,--zT~

V --xF

ilc

_

az

F

- 4 ; 06zlR,b

(14)

where RF includes the effect of adsorption on the porous matrix wall. q is the current of contaminant crossing the fracture/porous matrix boundary and may be shown to take the form q,-qQ$i.

ax

(15) x=b

where Q,is the porosity of the porous matrix and DM the molecular diffusion coefficient in the rock. In the porous rock matrix

ac,_h a+,

at-ET-

ilc

M ; blxlB

(16)

The boundary conditions are c&,0)

=0

C,(x,z,O)

=0

c,(O,t) = c,,

C,(b,z,O= C&t)

CF(DO,t)= 0

& ax

=o x=b

This set of equations is reminiscent of a reactor cell problem (albeit time-dependent) where the fracture is the fuel feeding particles into the rock matrix(moderator). This formulation provides useful information on the effect of fractures or fast legs in porous media. There are, however, some important questions regarding the formulation of eqn( 14) , the most important of which is how to distinguish the adsorbed fraction on the surface of the matrix from the actual solution CM in the matrix. The introduction of RF implies a surface layer of adsorbed atoms but the equation for CM includes no provision for the diffusion of these.

250

M. M. R. Williams

C Monte Carloji-acture generation

Some effort has been spent on the application of the Monte Carlo method to radionuclide and groundwater transport. Essentially, random numbers are used to generate a number of overlapping lines (fractures) according to a probability law specified by a mean fracture length and orientation distribution. Then, marked particles are tracked through these fractures according to some equation of motion determined by fluid flow and other processes such as adsorption. A large history of particle behaviour is accumulated and appropriate averages taken (de Marsily, 1986). This method is computationally time consuming and depends upon a knowledge of the underlying probability laws of fracture distribution and knowledge of the flow characteristics. If these are known accurately, the results are very good, if they are not, the uncertainty involved does not usually justify the effort needed. III

A TRANSPORT THEORY APPROACH

Let us look again at Fig 2, and imagine water flowing through the random matrix: then imagine a marked particle and follow it through all its twists and turns. It looks like a scattering process. There appear to be a series of rectilinear motions interspersed with sudden changes of direction. Looking more closely at an intersection, which we call a redistribution point, i.e. where two fractures intersect, the particle travels in a new direction as determined by the rock fracture orientation. This process of rectilinear motion followed by a change in direction is surely reminiscent of kinetic theory and suggests that one draws an analogy with neutron transport. The length of the rectilinear motion between redistribution points is a mean free path, and the probability of taking a new direction at a redistribution point is a scattering function. One can therefore define an angular concentration ~(r,~v,t)drdQdv which is the concentration of marked particles in volume element dr at r with speed between v and v+dv travelling along the rock fracture whose orientation is Q in the solid angle dQ. Using well known arguments one can derive a balance equation for C in the form (Williams, 1992)

1

R$+vR.V+vZ(R)+~ =

J J

dt2’ dv’v’C(CY)g(v

C(r,Qv,1)

+

v;Q’+ Q)C(r,iY,v',t)

+ S(r,R,v,t)

(17)

l/Z is the average distance between fracture intersections and may vary with orientation. g(...) dv dQ is the probability that a particle with speed v’travelling in the fracture whose orientation is Q’, will, at a distribution point, change speed to between v and v+dv and move into a fracture whose orientation is Q in dQ. To calculate Z(a) and g(...) requires a knowledge of the rock fracture angular distribution function, F(Q). Where F(d) dQ is the fraction of fractures whose orientation are Q in da. This can be obtained experimentally or by idealisations of the rock geometry. (There are some similarities here with the work of Ganapol (1989) on radiation transport through tree canopies). This transport equation has some clear advantages. It does not involve any explicitly random parameters; rather these are implicit in the form of F(R) and the equation is deterministic. Secondly it has the form of a time - dependent neutron transport equation for which many well-tried computer codes exist. Thus very little numerical development work is necessary. It remains to be seen whether the new equation tells us any more than the old one but it is likely to do so because it has more physical content and the parameters g(...) and Z(Q) are obtainable from geological data. We can rewrite the transport equation for a onedimensional problem as

Radionuclide

transport

251

where j.~=cos@8 being the angle that a fracture makes with the x-axis. Briefly we have assumed that the mean distance between fractures l/X does not depend on orientation, and no velocity change occurs when a particle passes a redistribution point, also the probability of a change of direction at an intersection depends only on the angle C2.iwhere i is the unit vector in the direction of the x-axis. These are not necessary assumptions but are useful to make. We have carried out two calculations to test the usefulness of the transport equation. In the first, an effective dispersion coefficient is calculated so that the classical solution matches as closely as possible the exact solution. This leads to a behaviour that is in close agreement with the experiments on scale dependence discussed earlier in section IIA. We note that the effective dispersion coefficient increases with distance from the source, in agreement with the experimental observations. We also note that it increases with time, thereby concurring with the results of the stochastic differential equation. This result indicates that the transport equation can explain the scale dependent effect and that it provides a viable alternative procedure. The second calculation is aimed at studying how the structure of the rock affects the concentration distribution. This is a problem that the classical equation is unable to deal with. In order to carry out the calculation, it is necessary to create a model of the rock. Thus we assume that the rock fractures are uniformly distributed as follows:

g(e)=&; -a181a

(19)

i.e. uniformly distributed between - a and +a. A small value of a means that the fractures are all orientated in much the same direction, whilst a large value means they are randomly orientated with a much greater spread of directions. It is easy to show that for the above distribution, I; = p7 (2~2-sin(%)) / 2a2, where p is the fracture density per unit area and is the average fracture length. We illustrate the results, by examining the fate of a source of radionuclides of constant height and duration 0.3 years for a range of values of a. Fig 5 shows the result. 0.01

0.008

0.006

..---.-*--

C(O.2,6.0)

- .w-..-.-.II.

C(O.4.6.0)

-----1-1-1

C(O.6,6.0)

-------

C(1.0.6.0)

-*-*-*-.-.

C(1.2,6.0)

-------

C(PI/2,6.0)

0.004

0.002

0 350

400

450

500

I

600

Distance from source in metres

Fig 5. Concentration as a function of fracture angle alpha

252

M. M. R. Williams

For the smaller values of a the fractures are nearly parallel and so the solute moves through the rock virtually by pure advection. As the orientation increases in randomness, so the distribution loses knowledge of the source and the rectangular nature tends to a smooth behaviour. Finally, when the fractures are randomly orientated over the half-sphere the distribution takes on a distinctly Gaussian shape. No other model of radionuclide transport, other than those based on Monte Carlo simulation, can predict such behaviour. Comparison of this model with the conventional advection-dispersion model has been made by Buckley et al (1993) and interesting and important differences are noted.

IV CONCLUSIONS

AND SUGGESTIONS

FOR FURTHER

WORK

We have described some of the basic physical problems and some of the current methods of modelling radionuclide transport. The limitations of those methods have also been stressed and, in particular, the inability of the advective-dispersion method to deal adequately with the scale dependent anomaly. The method of stochastic differential equations goes some way to giving an explanation of this phenomenon but is impractical for large scale calculations. The double porosity method is described and answers some of the problems raised by the presence of fractures, but its underlying derivation is obscure. A new method is proposed, based upon an analogy with neutron transport theory, in which fractures are analogous to neutron mean free paths and intersections of fractures analogous to scattering. This method explains the scale dependence, with no assumptions, and also has the marked advantage of not requiring much numerical development, because there are several excellent computer codes available for time dependent neutron transport, which can easily be adapted to the radioactive waste problem. Finally, we note a number of areas in which further effort is needed. 1. Interaction of the diffusion and sorption processes in fractured rock, e.g. does matrix diffusion retard and at what rate does the system reach equilibrium ? 2. The validity of reversible sorption. Is it legitimate to assume that adsorption and desorption follow the same isotherm or is there some effect of hysteresis ? 3. Dynamics of groundwater with adsorption and desorption, e.g. the effect of a moving saline front leading to spatial changes in I&, 4. Coupled flow and transport problems. Currently the groundwater flow is obtained independently of the nature of the solute. There may be interactions, especially if heat is generated. 5. The way in which the rock matrix is deformed by the fluid moving through it. This involves a coupled stress, porosity and flow problem.

ACKNOWLEDGEMENTS The author is grateful to Dr M.C. Thome for much valuable advice. He also wishes to thank Professor Ulf Lindblom for permission to reproduce Fig 1 and to Professor Ivar Neretnieks for permission to reproduce Fig 2. Acknowledgement is also given to Atomic Energy of Canada Ltd from which Fig 3 was obtained.

Radionuclide transport

253

REFERENCES Brenner, H. (199 1) “Macrotransport Processes”, Eighteen years of colloid and surface chemistry. The Kendall Address, 1973- 1990, Edited by T. Fort and K.J. Mysels. Buckley, R.L., S.K. Loyalka and M.M.R. Williams, (1993) “Numerical studies of transport in porous media III: radioisotope migration in fracture angle formations”, Ann. Nuclear Energy 20, 569. de Marsily, G.(1986) Quantitative HydrogeoZogy , Academic Press. Ganapol, B. (1989) “Radiative transfer in dense plant canopies with azimuthal symmetries”, TTSP 18, 475. Sudicky,E.A. and E.O. Frind,( 1982) “Contaminant transport in fractured porous media: analytical solutions for a system of parallel fractures”, Water Resources Research 18, 1634 Williams, M.M.R.( 1992) “A new model for describing the transport of radionuclide through fractured rock”, Ann. Nuclear Energy 19, 791. Further Reading not explicitly mentioned in the text:

J. Bear, Dynamics of Fluids in Porous Media, Dover Publications (1988) J. Bear, Hydraulics of Groundwater, McGraw Hill Book Co. (1979) A. Verruijt, Groundwater Flow, MacMillan Press (1982) K.B. Krauskopf, Radioactive Waste Disposal and Geology, Chapman and Hall (1988) U. Lindblom and P. Gnirk, Nuclear Waste Disposal, Can we rely on bedrock ?, Pergamon Press (1982) M.M.R. Williams, “Stochastic Problems in the transport of radioactive nuclides in fractured rock”, Nucl. Sci. and Engng. 112, 215 (1992) Professor Nils Sjiistrand: I owe a personal debt of gratitude to Nils Sjiistrand. Many years ago as a young research student working on my Ph D thesis, I was inspired by his paper in Arkiv Fysik on the variation of extrapolation distance with buckling using one-speed transport theory. That work guided me and suggested extensions to energy dependence that were crucial in gaining a better understanding of the pulsed neutron experiment. Thank you Nils.