Direct simulation of heterogeneous diffusion and inversion procedure applied to an out-diffusion experiment. Test case of Palmottu granite

Journal of Contaminant Hydrology 93 (2007) 21 – 37 www.elsevier.com/locate/jconhyd

Direct simulation of heterogeneous diffusion and inversion procedure applied to an out-diffusion experiment. Test case of Palmottu granite P. Sardini a,⁎, J.C. Robinet a,b , M. Siitari-Kauppi c , F. Delay a , K.H. Hellmuth d a

c

HYDRASA Laboratory (Hydrogeology, Clays, Soils and Alterations), UMR 6532 CNRS, University of Poitiers, 40 avenue du Recteur Pineau, 86022 Poitiers Cedex, France b ANDRA, 1/7 rue Jean Monnet, Parc de la Croix Blanche, 92298 Châtenay-Malabry cedex, France Radiochemistry Laboratory, Department of Chemistry, University of Helsinki, P.O. Box 55, FI-00014 Helsinki, Finland d Finnish Centre for Radiation and Nuclear Safety (STUK), P.O. Box 14, FI-00881 Helsinki, Finland Received 3 February 2006; received in revised form 13 November 2006; accepted 15 January 2007 Available online 26 January 2007

Abstract An out-diffusion laboratory experiment using a non-reactive tracer was fitted using the Time Domain Diffusion (TDD) method. This rapid particle tracking method allows simulation of the heterogeneous diffusion based on pore-scale images and local values of diffusivities. The superimposed porosity and mineral 2D maps act as computation grids to condition diffusion pathways. We focused on a Palmottu granite sample, in which the connected pore space has a composite microstructure with cracks linking microporous minerals and is above the percolation threshold. Three main results were achieved: (i) When compared to the fitting obtained with one coefficient (best mean square residual R = 1.6 × 10− 2), diffusion is shown to be suitably characterised with two coefficients related to cracks and microporous minerals (best R = 6.5 × 10− 4), (ii) rather than imposing a local apparent diffusion coefficient Da independent of the local porosity Φ, a best fit is obtained by applying Archie's relationship Da = D0 × G with G = Φm to each pixel of the calculation grids (G is the geometry factor, D0 is the diffusion coefficient in free fluid, and m is Archie's exponent), and (iii) the order of magnitude of the fitted diffusion coefficient or Archie's exponents (m = 0 for microcracks and m = 1.82 for microporous minerals) is physically realistic. © 2007 Elsevier B.V. All rights reserved. Keywords: Matrix diffusion; Out-leaching experiments; Archie's law; Granite; Porosity; Inversion modelling; Finland; Nuclear waste disposal

1. Introduction In the context of an underground nuclear waste disposal, the possible radionuclide migration through crystalline rocks has been studied extensively from in situ or laboratory experiments. A given volume of crystalline rock in the field is structurally characterised by its

⁎ Corresponding author. Tel.: +33 549453828; fax: +33 549454241. E-mail address: [email protected] (P. Sardini). 0169-7722/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2007.01.011

fractures and matrix. Fractures represent complex waterconductive structures, spanning scales from a few centimetres to several decametres (Mazurek et al., 2003). The matrix refers to zones of intact or altered and/or deformed rock between the fractures. Two domains can be observed in the matrix: (i) Regions located near some fractures, characterised by non-stationary spatial characteristics due to deformation gradients (fault gouge, cataclasite, mylonitic bands…) and/or alteration gradients;

22

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

(ii) Spatially stationary media, often named “fresh” rock that may also enclose large pervasively altered volumes connected to conductive fractures. Diffusion within the rock matrix is considered to be one of the predominant processes for the migration of radionuclides (Bear et al., 1993). Diffusion in lowly permeable granitic rock is an important issue in the disposal of spent fuel from nuclear powerplant, since several countries are planning to locate the repositories in granitic formations. Matrix diffusion represents one of the main retardation factors affecting migration of pollutants (Neretnieks, 1980; Carrera et al., 1998). At present, our knowledge of the complex pore structure in the matrix is poor, making it difficult to analyse accurately the effect of diffusion. Thus, laboratory diffusion experiments or field experiments in heterogeneous crystalline matrix are somewhat quantitatively hypothetical (Skagius and Neretniecks, 1986; Jacob et al., 2003; Vilks et al., 2003). Recent studies have modelled heterogeneous matrix diffusion from statistical distributions of diffusion coefficients in the rock, which control the residence time distribution in the rock matrix (Haggerty et al., 2000; Tidwell et al., 2000; Vilks et al., 2003). It was suggested by Haggerty (2002) that the distributions of diffusion coefficients should be conditioned on pore space connectivity. Hellmuth et al. (1993) have developed an impregnation technique suited to the structural analysis of low porosity crystalline matrices. Autoradiographic imaging of slices of impregnated rocks allows us to visualise and to quantify the rock connected porosity at the scale of a core (Siitari-Kauppi, 2002). Each pixel of a 2D image is characterised by a known porosity, making it possible to calculate porosity distributions used to analyse diffusion experiments (Johansson et al., 1998). The aim of this paper is to point out the influence of pore connectivity on diffusion in heterogeneous crystalline matrices. We demonstrate how to characterise the distribution of apparent diffusion coefficients inside the matrix from an adapted simulation procedure. A cylindrical core of the Palmottu granite, Finland, was submitted to an out-leaching experiment using a non-reactive tracer (14C-methylmethacrylate, 14C-MMA). The same sample was then saturated with 14C-MMA, and the porosity maps of eight sections were obtained from autoradiographs. The out-leaching process was modelled over porosity maps from the autoradiographs using the Lagrangian method of Time Domain Diffusion (TDD) for the solute transport (Delay et al., 2002; Sardini et al., 2003). An inversion procedure was applied to automatically adjust the experimental out-diffusion curve by varying diffusion coefficients on subareas delineated from the spatial distribution of

porosity and primary minerals. To optimise the computation time the Jacobian matrix used in the inversion process was derived analytically, as proposed by Delay and Porel (2003). The method was slightly modified to account for the apparent diffusion coefficient that obeys a power law relationship with the local porosity. Particular attention was paid to the image conditioning before inversion. Images for inversion were thresholded in porosity classes, and we have verified the minimal number of classes needed to avoid numerical degradation of the porous network. For computational optimisation, we also needed to seek the most appropriate procedure for thresholding the conductive pixels, and more importantly, for determining the optimal pixel size. 2. Methodology 2.1. Out-diffusion experiment, porosity and mineral maps The sample studied in this work (sample reference R388-33m) is a core of 45.5 mm in diameter and 50 mm in length from the Palmottu granite, Eastern Finland. It is a coarse-grained porphyritic granite containing mainly centimetric alkali feldspar, plagioclase, quartz and biotite (Siitari-Kauppi et al., 2003; Oila et al., 2005). The average porosity of the whole rock matrix is 0.5 ± 0.1% regardless of the measurement method, 14C-MMA or water saturation. The core was first dried for 1 month at 80 ± 5 °C. It was then saturated with 14C-MMA resin for 3 weeks under vacuum (10 Pa). 14C-MMA is a geochemically nonreactive tracer. It was employed to avoid the difficulties inherent in the distinction between pure matrix diffusion process and sorption processes (Carrera et al., 1998). Starting from a constant concentration of 14C in the connected pore space, an out-diffusion experiment was carried out on the same sample. Radioactive MMA was leached from the core in 400 ml of inactive MMA for 3 months. Samples (0.5 ml) from the external solution were collected regularly and analysed by scintillation counting. The activity of the total external volume versus the square root of time is plotted in Fig. 1. This experimental leaching curve (Fig. 1) is slightly different from the initial curve provided by Sardini et al. (2003); some points have been removed from the initial curve in order to smooth it. It is assumed that using the initial curve as an experimental input would not modify the inversion results that are presented below. The same sample was saturated with 14C-MMA for 3 weeks for the second time. After resin polymerisation using a 60Co source, the core was sawed normal to its long axis to get eight flat cylindrical slices. In brief, the

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

23

Fig. 1. Experimental leaching curve of R388-33m Palmottu granite cylindrical sample. The sample dimensions are 45.5 mm in diameter and 50 mm in height.

autoradiographic technique provides 2D pictures of the 3D connected porosity network (Hellmuth et al., 1993; Siitari-Kauppi, 2002). Examples of three autoradiographs from these slices are shown in Fig. 2. All of them display the same porous sets, i.e. a composite microstructure with cracks of micrometric aperture connecting porous patches (Oila et al., 2005). Fig. 2 also shows the total 2D connectivity of this pore network. This means that this network forms an “infinite cluster” (Stauffer and Aharony, 1994) clearly above the percolation threshold, and that the diffusion process can reach any point in this network. Cracks are visible at the grain boundaries around alkali feldspar and quartz. Opened cleavage planes of alkali feldspar are also very visible. The porous patches are mainly hematised plagioclase and chloritised biotite. The spatial distributions of porosity slightly differ between autoradiographs, due to variations in mineral contents between the slices. These variations result from the poor

spatial representativity of the sample because of the presence of centimetric alkali feldspar phenocrysts. Slice 2 (average porosity Φ = 0.47%) shows an “average” porosity pattern, whereas slice 5 (Φ = 0.38%) contains mainly microcracks because it crosscuts a large alkali feldspar phenocryst. Slice 7 (Φ = 0.49%) contains many porous hematised plagioclase crystals which are well connected to the border of the sample. In the following investigations, we will specifically study these three slices because of these differences in terms of pore network connectivity. A staining technique was applied to these three slices to determine the spatial distribution of primary minerals. This technique is based on chemical etching and staining of the main primary minerals of the rock (Muller, 1967). It proved particularly efficient for revealing the mineral distribution of coarse-grained Palmottu granite. Mineral maps were obtained from the stained surfaces by applying simple colour-thresholding methods (Sardini et al., 1999).

Fig. 2. Autoradiographs from three slices of the sample R388-33m (exposure time 5 days, diameter of the slices 45.5 mm). These images were the basis for diffusion simulations and inversion of the observations shown in Fig. 1.

24

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

Each colour corresponds to a single primary mineral specie (Fig. 3a). Mineral and porosity maps were numerically superimposed using PhotoShop™ software. They constitute the two layers used in the quantification of heterogeneous diffusion in Palmottu rock. 2.2. Time domain diffusion (TDD) simulation TDD was developed specifically for the fast calculation of diffusion processes over digitised autoradiographs (Delay et al., 2002; Sardini et al., 2003). Digitised autoradiographs constitute the numerical computation grids for diffusion calculations. For a given pixel i, we consider that the porosity Φi is known from the autoradiograph, and the apparent diffusion coefficient Di is the sought parameter. The diffusion equation in heterogeneous and isotropic media can be written as: Að
CÞ ¼ jd ð
Dd jCÞ At

ð1Þ

where Φ is the porosity [−], D the apparent diffusion coefficient [L2T− 1] and C the concentration [ML− 3]. Using a finite volume formalism, Eq. (1) may be rewritten as: Cj −Ci ACi X
i Vi ¼ Sij ð
DÞij At Lij j

porosity and apparent diffusion coefficient between pixels i and j. A partial derivative with respect to time is a material derivative in the Lagrangian context. Thus, Eq. (2) can be rewritten in the form: " # X X dCi ¼− bij Ci þ bij Cj ð3Þ
i Vi dt j j where bij =Sij(ΦD)ij /Lij. For a regular network of pixels of size 1 on an edge, Sij =12, Lij =1 and bij [L3T− 1] =1(ΦD)ij. Eq. (3) is equivalent to:

/i Vi P

bij

dðlogðCi ÞÞ ¼ −dt þ

j

X Cj bij dt Ci j

P

bij

ð4Þ

j

The second term on the right-hand-side shows that expression (4) may be separated into j transition equations P with an occurrence probability Pij ¼ bij = j bij. Thus, the diffusion time from i to j can be calculated as a value drawn from an exponential distribution and simulated by: Ui Vi tiYj ¼ − P logðu01 Þ bij

ð5Þ

j

ð2Þ

indexes i, j refer to pixel i having four pixels j (or six voxels j in 3D) which are direct neighbours of i. Vi is the volume of pixel i, Cλ the mean concentration in pixel λ (λ =i, j), Lij the lag distance between the centres of pixels i and j, Sij the total area of the interface between i and j, and (ΦD)ij the harmonic mean of the product of

u01 being a random number from a uniform distribution between 0 and 1. A model able to mimic any diffusion experiment is obtained duplicating the algorithm N times, N being a sufficiently large number of particles to attain statistical stability of the calculated concentrations. This TDD method is much more rapid than the classical random walk method because the particles are forced to move a distance l at each jump. This is not the case with

Fig. 3. Example of slice 2: the mineral map (a) is superimposed over the simplified autoradiograph (b) with a resolution of 200 dpi (QZ = quartz, PL = plagioclase, KF = alkali feldspar and dark spots are biotite crystals). The map of the porous subsets (c) with fissures and porous patches were deduced from images (a) and (b). Black corresponds to the free fluid around the sample, white corresponds to impervious pixels, dark grey to porous patches, and light grey to fissures. Diameter of the slice 45.5 mm.

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

the classical random walk method, because the particles have to undergo numerous jumps within the same pixel for the method to be accurate. 2.3. Inversion of the diffusion process Inversion of direct diffusion modelling was undertaken in order to identify the apparent diffusion coefficient Da of each pixel in the system. Inversion procedures are useful because direct calculations of diffusion are protracted and also because foreseeing the macroscopic diffusion coefficient of any natural material is not straightforward. We have used the inversion method developed by Delay and Porel (2003) to determine automatically the apparent diffusion coefficients for different porous subsets. The method is based on the minimisation of a least square objective function using a Gauss–Newton algorithm, the terms of the Jacobian matrix being analytically derived. Consider a vector d of size n representing the set of diffusion coefficients to be determined, and a vector ε of size p related to the sampled errors between observation and simulation. The closeness criterion between observation and simulation at iteration k, i.e. the objective function, is written: F k ¼ 1=2ðek ÞT ek

ð6Þ

with T the transposition operator. The Taylor development up to the second order of F at iteration k + 1 and with respect to d is: F kþ1 ¼ F k þ jF k d ðd kþ1 −d k ÞT 1 þ ðd kþ1 −d k ÞT d j2 F k ðd kþ1 −d k Þ 2

ð7Þ

∇F and ∇2F are respectively the gradient vector (∂F / ∂di) and the hessian matrix (∂2 F/(∂di∂dj)). Minimisation of F k+1 is reached for ∇F k+1 = 0. Thus, calculating the derivative of F k+1 in Eq. (7) with respect to d k+1 and cancelling it yields an expression for calculating d k+1 : jF k ¼ −j2 F k d ðd kþ1 −d k Þ

ð8Þ

∇F and ∇2F depend on Jε, the p × n Jacobian matrix of the errors. Each term in Jε is defined by the partial derivative of the errors εα with respect to the parameter dβ (α : 1..p and β : 1..n) : ∂εα / ∂dβ. The Gauss–Newton procedure uses the following exact expression for the gradient vector and a first order approximation for ∇2F: jF ¼ JeTk d ek ; j2 F ¼ JeTk d Jek

ð9Þ

The terms of Jε can be approximated by the perturbation method; however this technique requires n direct

25

calculations of the diffusion problem at each iteration (Delay et al., 2002). Delay and Porel (2003) proposed a faster method based on the analytical calculation of Jε. The method requires only one direct calculation at each step to compute all the terms of the Jacobian matrix. The terms ∂εα / ∂dβ of Jε can be rewritten as the product of two partial derivatives (∂εα / ∂t) × (∂t / ∂dβ). The first term is approximated by a finite difference near the sampling point α: Aea =Atc

eaþ1 −ea taþ1 −ta

ð10Þ

The term ∂t / ∂dβ is determined by analytical calculation. If one considers a particle making an elementary jump from pixel i to pixel j, the transition time ti →j (5) is a function of l, Dλ and Φλ with λ = i,j. Thus for each particle, where parameter d refers to diffusion coefficient D, ∂ti →j / ∂dβ can be expressed as: ð11Þ 8kaf j ¼ neigh: of ig; ðAtiYk =Adb Þ=lnðu01 Þ " # 2 2 X ðUDÞij Ui 1 ¼ P 2 dðDi −db Þ D2i Ui j 2 ðUDÞij j " #! X ðUDÞ2ik þ dðDk −db Þ 2 Dk Uk kaj δ() is the Kronecker function (δ(0) = 1; δ(x ≠ 0) = 0). Calculation of the transition time derivative does not require generating again the random number u01 used in calculating the transition time. This means that the transition time and all its associated derivatives with respect to dβ (β : 1..n) can be calculated simultaneously. This accelerates the calculation considerably. For a given particle, the transition time derivatives are calculated at each jump, and are cumulated in order to find the time derivative P of the complete particle walk: Atparticle = Ad ¼ i;k AtiY =Ad. Finally, ∂t / ∂d is obtained merely by averaging all the time derivatives of all the particles leached at the sampling time α: 1=Nparta  P particle ðAtparticle =AdÞ. If a given particle has not leached out of the sample by sampling time α, it does not contribute to the measured concentration, and does not have to be included in the calculation of b∂t / ∂dβ N. 2.4. Conditioning diffusion coefficients from rock microstructure The above inversion procedure was applied to outdiffusion experiments performed on Palmottu granite. Digitised autoradiographs were made of hundreds of thousands of pixels. To minimise the number of sought

26

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

diffusion coefficients, either (i) a unique diffusion coefficient or (ii) two diffusion coefficients were used for all conductive pixels, one for all the pixels of the porous patches and one for the fissures. In the later case, the porous patches were defined by all pixels associated with hematised plagioclase and chloritised biotite from the superimposed mineral map. All other pixels were assigned as cracks: grain boundary microcracks, quartz intragranular fissures or open cleavage planes in K-feldspar phenocrysts. In addition, another way to constrain the system was tested, this time accounting for the local dependence of Da on porosity. This dependence is often referred to as the “Archie's law”, which links the apparent diffusion coefficient Da to the porosity by a power law: Da ¼ GD0 ¼ D0 U

m

ð12Þ

where G represents the geometry factor, D0 the diffusion coefficient in free fluid and m a power exponent obtained experimentally. This empirical expression has often been discussed in the literature (Worthington, 1993). Sato (1999) advocated the worth of Eq. (12) on the basis of experiments on igneous rocks, while Parkhomenko (1967) proposed a slightly different relationship between G and Φ: G = 0.71×Φ0.58. Other authors considered that the effects of pore space heterogeneity are not well reproduced by this macroscopic relationship (Skagius and Neretniecks, 1986; Siitari-Kauppi, 2002). We used Eq. (12) in this study at the pixel scale, and we determined the best power law exponents by the inversion method depicted above. The terms of the Jacobian matrix become ∂εα / ∂mβ, and the transition time derivative ∂ti →j / ∂mβ is linked to ∂ti →j / ∂dβ (cf. expression (11)) by: AtiYj =Amb ¼ Di lnðUi ÞAtiYj =Adb

ð13Þ

It is clear from Eq. (13) that introducing a power law linking Da to Φ does not modify the general inversion procedure. The same number of power law exponents has been used to adjust the experimental curve, i.e. one or two according to the sorting of pervious pixels in patches and cracks. 2.5. Out-diffusion modelling from 2D to 3D In order to inverse the 3D experimental leaching curve using 2D images of the pore network, simulated curves were normalised. First, reservoir concentrations C(t) were normalised in relation to the final concentration Cf reached at equilibrium. Final concentrations from experiment and simulation were determined according to a mass-balance calculation using (i) the initial activity of 14C-MMA in the pore space, (ii) the reservoir and sample volumes, and (iii)

the connected porosity of the sample. Those concentrations obtained by simulation are represented by a number of particles, whereas those obtained experimentally are given in Bq/ml. These activities were converted into a number of particles by setting the experimental equilibrium concentration to the final concentration obtained by simulation. Second, square roots of times were normalised against volume-to-surface ratios (V/S) of both 2D images and 3D samples. For 2D images, the V/S ratio is equal to the surface-to-perimeter ratio. The classical “square-root relationship” (14) provided by Hespe (1971) relates the slope of the curve (normalised concentration − square root of time) to the apparent diffusion coefficient Da: pffiffiffiffiffiffiffiffiffiffiffi pffi CðtÞ V  ¼ 2 Da =p  t Cf S

ð14Þ

This relationship, written for both 2D and 3D diffusion, was used to convert the square root of times from 2D to 3D according to both V/S ratios (15): pffi pffi ðV =SÞ3D t j3D ¼ t j2D  ðV =SÞ2D

ð15Þ

Expression (15) allows a direct comparison between two leaching experiments performed with two samples of the same rock having different external geometries. After that step, 2D and 3D curves were directly compared by the inversion algorithm. The least square residuals R obtained from experimental (3D) and P simulated 2 (2D) leaching curves were given by R ¼ a¼1 to p ea = P 2 a¼1 to p Ca ðexpÞ, where Cα(exp) represents the experimentally (3D) measured concentrations of the p observation points. A fitting was considered as reliable if R was close to 10− 3. In order to validate the proposed approach, two simulation tests related to virtual homogeneous and heterogeneous media are presented. First, an out-leaching simulation was conducted on a cubic and homogeneous rock sample (edge size 2.2 mm, Da = 10− 11 m2 s− 1, Φ = 1%) (Fig. 4). Second, a leaching curve related to a 2D square section (2.2 × 2.2 mm2) of this cube was calculated without any normalisation of the square root of times. Third, this last curve was normalised according to Eq. (15). It was observed that this normalised curve was well superimposed over the curve obtained from the cube, even though a slight shift between them was visible (see Fig. 4). Finally, inversion of the curve calculated from the cube using a 2D square section indicated an adjusted value of Da = 9 × 10− 12 m2 s− 1 (R = 8.6 × 10− 4). This value is very similar to Da = 10− 11 m2 s− 1, proving the reliability of the fitting methodology for the homogeneous case.

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

27

Fig. 4. Out-leaching curve modelled for a homogeneous cubic sample (Φ = 1%, Da = 10− 11 m2 s− 1). Fitting of this 3D curve was achieved with the help of the proposed inversion procedure using a square section of the sample for diffusion modelling.

Second, a heterogeneous cubic sample was investigated. It was built as a regular and periodic assemblage of small cubic porous patches linked by fracture planes (see Fig. 5 for a 3D view of sample geometry, and see Table 1 giving the parameters of fracture and porous patches). This geometry tries to mimic the composite microstructure with cracks connecting porous patches observed on the autoradiography of the Palmottu granite. Any slice of this sample gives a connected 2D pore network. Voxels which were neither fractures nor patches were considered as impervious (in white in Fig. 5).

Initially, three out-leaching simulations were successively performed on this volume using three Da for the patches, which were chosen as 5 × 10− 10, 5 × 10− 11 and 5 × 10− 12 m2 s− 1 (Table 2). An example of modelled outdiffusion curve using the 3D geometry is presented in Fig. 7; in this graph, this out-diffusion curve was normalised to a 2D square slice (1 cm2) according to the inverse of Eq. (15). The curve obtained using the volume is comparable to a real 3D experiment. In order to inverse diffusion curves originated from 3D geometry, four 2D square slices named A, B, C, and D

Fig. 5. 3D perspective of the virtual heterogeneous sample used to validate the diffusion and inversion method. Parameters of the different porous subsets (patches, fractures) are given in Table 1. White corresponds to impervious voxels.

28

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

Table 1 Parameters of the virtual sample (Fig. 5) used for heterogeneous diffusion modelling Sample parameters

Patches parameters

Fractures parameters

Size 1003 voxels Voxel size 0.1 mm

Size 133 voxels Da (m2 s− 1) 5 × 10− 10 5 × 10− 11 5 × 10− 12 Porosity 10%

Thickness 1 voxel Da (m2 s− 1) 10− 9

Porosity 1%

were tested (see Table 2 and Fig. 6). Similarly to the autoradiograph and mineral map, one slice characterises the spatial distribution of pores and minerals. These four slices try to represent possible variations of the pore network for random slices of this periodic volume. They were built manually, and do not represent real 2D slices of the volumic sample. We emphasise that studying real random slices of this periodic volume would have been more complex than studying random slices of isotropic random porous media, which would have been statistically similar. Thus, the results presented below are only indicative of some constraints when studying diffusion within heterogeneous media, and are quite difficult to transpose to real 3D porous media. The porosity of a slice and the connectivity of porous patches with the border of a slice were selected as the two variables influencing diffusion (Table 2). The connectivity of porous patches with the border of a slice was represented by Pe/Pt ratio, Pe representing the patch volume that is connected to the border of the slice and Pt being the total patch volume. Sardini et al. (2003) suggest that this ratio controls the position of an inflexion point in an out-leaching curve related to such composite heterogeneous media. This point corresponds to the end of the leaching of the patches connected to the border of the slice. Compared to the parameters of the tested volume, slice A has the same porosity (1.37%)

and the same Pe/Pt ratio (0.58), slice B has a slightly lower porosity (1.08%), slice C has a higher porosity (1.96%), and slice D has a lower Pe/Pt ratio (0.41). Furthermore, the fracture density of these four slices was chosen equal to the density viewed in a section parallel of a face of the cube. Note that an important difference between 2D and 3D representations is that fractures of the studied slices in 2D are not directly connected to the border of the slice, whereas in 3D, fractures are connected to the sample's border. Examples of inversion of a 3D out-diffusion curve are presented in Fig. 7 (Da = 5 × 10− 11 m2 s− 1 for the porous patches). The sought parameters were the apparent diffusion coefficients of patches and fractures. Concerning the fractures, all inversion tests yielded the same Da as the 3D one (10− 9 m2 s− 1). The diffusion coefficients of patches and related least square residuals obtained by inversion are given in Table 2. These results demonstrate that porosity, and more importantly, patch connectivity with the border of slice influence the inversion procedure. Indeed, using slices A and B provides a satisfying fitting of the 3D curves (Fig. 7, Table 2), the fitted diffusion coefficients being close to the 3D ones (Table 2). By contrast, it was never possible to get a good fitting using slice D, because in terms of normalised concentration, the inflexion point of the leaching curve was located too low, in the corresponding Pc/Pt ratio of slice D (0.41). Moreover, using slice D, the diffusion coefficients resulting from inversion were far from the 3D ones. Fitting using slice C was better than fitting using slice D, even if it was not acceptable in terms of least square residuals. However, the found parameters using slice C were close to the 3D ones. To conclude, the out-diffusion process in a slice of a heterogeneous medium can be relatively difficult to handle owing to its dependence on the connectivity of porous subsets with the border of the slice. Thus, for the studied geometry, the slice chosen to represent the 3D pore space must have nearly the same Pc/Pt ratio than the 3D sample. Furthermore, porosity of the tested slice

Table 2 Patch connectivity and porosity of the four slices presented in Fig. 6, and results of inversions for these four slices and for different parameterisation of diffusion

Porosity Pc/Pt Da (patches) (m2 s−1)

Volume

Slice A

Slice B

Slice C

Slice D

1.37% 0.58 5 × 10− 10

1.36% 0.58 2.3 × 10− 10 R = 1.9 × 10− 3 3.6 × 10− 11 R = 1.3 × 10− 3 7 × 10− 12 R = 3.1 × 10− 3

1.08% 0.58 2 × 10− 10 R = 1.3 × 10− 3 3.5 × 10− 11 R = 9.5 × 10− 4 5.1 × 10− 12 R = 1.1 × 10− 3

1.96% 0.58 4.9×10− 10 R = 9.6 × 10− 3 4.4 × 10−11 R = 5.7 × 10− 3 4.4 × 10− 12 R = 7.7 × 10− 3

1.38% 0.41 −

5 × 10− 11 5 × 10− 12

10− 10 R = 1.1 × 10− 1 1.5 × 10− 11 R = 4.4 × 10− 3

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

29

Fig. 6. These four bidimensional representations of the 3D sample presented in Fig. 5 were constructed in order to test the effect on the inversion procedure of both porosity and patch connectivity with the border of the slice. For each slice, patch connectivity (Pc/Pt) and porosity are provided in Table 2.

must be close to the volumic porosity. When using an adequate slice, it is demonstrated that the employed inversion approach allows obtaining a good approximation of the local diffusion coefficients related to each porous subset of the 3D sample. When adjustment was satisfying, the deviation between apparent diffusion coefficients related to the volume and those determined by inversion using a 2D slice was never greater than a factor 2.5 (Table 2). 3. Application to samples from Palmottu granite 3.1. Image conditioning Even with the TDD method, direct simulation of diffusion over a large image is a cumbersome task. Moreover, often more than twenty iterations are required to reach convergence in an inversion procedure. Thus, image pre-conditioning is applied to decrease the computation time.

We first investigated the sensitivity of the diffusion calculation to the number of porosity levels in the porosity maps. The number of porosity levels did not influence the calculation time. Nevertheless, a test was performed to assess the acceptable degree of simplification of the maps. Direct calculations of out-diffusion were computed over an image of a slice of Palmottu granite, with the number of porosity levels ranging from 3 to 35 (Fig. 8). Above 12 porosity levels, the out-leaching curves were very similar, and were closely conformed to those at 8–12 levels. This means that with reasonable simplification, i.e. 8–16 levels, one gets a good approximate image of the porosity distribution. Note, however, that this result may be theoretically valid only for the porous network of Palmottu granite, and cannot be generalised because the statistical distribution of porosities within an image influences the calculation stability. In subsequent work, 16-level porosity maps were used. The thresholding of a 2D connected porous network must be undertaken carefully because it defines the 2D

Fig. 7. Example of an out-leaching curve modelled for the volume presented in Fig. 5, and computed with Da = 5 × 10− 11 m2 s− 1 for the patches (square dots). Solid lines represent simulated best fit curves using 2D slices A, B, C and D. The resulting fitting parameters and related least square residuals are given in Table 2.

30

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

Fig. 8. Out-leaching curves modelled for the same image of a slice, but with a number of porosity levels ranging from 3 to 35. The leaching curves show stability above twelve porosity levels.

connectivity of the network. Using the 14C-PMMA method, impervious pixels are defined from a porosity level also called the “background” level. This limit is measured on a zone of the autoradiographic film considered as unexposed (Siitari-Kauppi, 2002), and is often selected as the maximal grey level found on the scanned image. However, the use of background as a single threshold may yield underthresholding of impervious pixels. Indeed, some impervious zones in the autoradiographic film were darkened by the lateral electronic emission originating from neighbouring pores and present slightly lower grey levels than the background. So, we employed the watershed method based on the gradient image calculation and on flooding algorithms (Coster and Chermant, 1989). Routines from the Micromorph™ software were used as follows. The autoradiographs were initially filtered by an elementary opening, and gradient images determined from these filtered images. On each image, marker zones were needed to initiate the flooding process affecting the gradient image. These zones were chosen as the union of extended maxima of the filtered autoradiograph and of rough thresholding of the porous patches. The final segmentation was therefore influenced by two variables which were adjusted manually: (i) the order of the extended maxima, and (ii) the threshold level used to extract the porous patches. This method provided satisfying results, with each conductive pixel being connected to at least one conductive neighbour by a border (e.g. see Figs. 2a and 3b). Thresholding minimised the number of conductive pixels, thus optimising the computing time. The most important step of image conditioning was to determine the optimal pixel size, l. Pixel size has a radical

influence on computation time. If 1 is low then the average transition time bti →jN is small, and the number of transitions needed to travel a given distance is high. From Eq. (5), ti →j is a quadratic function of the pixel size. Thus, the computation time for a given number of injected particles will be approximately a function of 1 −2. Imposing a weak definition on the porosity maps (high value of 1) significantly reduces the computation time. However, the use of weak definition images also yields some changes in the connectivity of the porous network. Conversely, the use of a small pixel size preserves the 2D connectivity of the network, but the computation time becomes too long to solve the inverse problem. The effects of image definition on 2D connectivity and simulation of diffusion were investigated for 1 ranging from 63.5 μm (400 dpi) to 254 μm (100 dpi). All images were built from the 400 dpi-scanned autoradiographs by reducing the image definition using the bicubic algorithm. Two connectivity indexes C1 and C2 were calculated for the digitised network of conductive pixels: C1 ¼

Nc Nc −Nde ; C2 ¼ Ntot Nc

ð16Þ

Nc is the number of conductive pixels, Ntot is the total number of pixels composing the sample (conductive + impervious), and Nde is the number of pixels detected within 2D dangling-ends. The dangling-end pixels controlled by 1 or 10 pixel-sized necks were successively detected (Nde1 and Nde10) using the algorithm by Sardini et al. (1997). Fig. 9 presents the evolution of C1 and C2 (computed for Nde1 and Nde10) as a function of image definition. C1 remains almost constant from 400 to 200 dpi, but increases clearly from 200 to 100 dpi. For

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

31

Fig. 9. Evolution of the 2D connectivity with respect to an image resolution from 100 to 400 dpi. C1 increases abruptly above a 200 dpi resolution. C2 varies weakly within the same resolution range.

the low definitions, the pixel size is larger than the mean aperture of the fissures, which yields an overestimate of the ratio of conductive pixels to the total number of pixels. C2 increases continuously with increasing pixel size. However, these variations are small because the number of dangling-end pixels is small when compared to the total number of conductive pixels. C2 is almost stable from 400 to 200 dpi when calculated with Nde1 and varies from 0.91 to 0.93 when calculated with Nde10. The effects of image definition were also assessed on numerical simulations of out-diffusion (Fig. 10). For these tests, unique values of the local diffusion coefficients were assigned to the conductive pixels. The leaching curves are very similar in shape, but some slight variations are conspicuous. These differences are linked globally to image definition even if slight discrepancies are visible. For example, the two similar upper curves in Fig. 10 correspond to 200 and 400 dpi resolution images, although globally the lower the definition the lower the diffusion velocity.

Considering the results given by the geometrical 2D connectivity and the simulated leaching curves, a 200 dpi resolution (l = 127 μm) was selected for the inversion modelling. This resolution appears to be a good compromise between reasonable computing time and good connectivity preservation of the porous network. 3.2. Inversion of out-diffusion experiment Out-diffusion experiment on sample R388-33m was inverted by 2D simulations on slices 2, 5 and 7 (Figs. 2 and 3). The inversion process yielded a valuable solution within 5 to 30 iterations. The sought parameters and related least square residuals are listed in Table 3. The following comments may be drawn from the results in Table 3 and in Fig. 11: • A single variable, i.e. just one Da or one m for all conductive pixels, is never enough to fit an experimental curve whichever slice is considered. Fitting

Fig. 10. Out-leaching curves as a function of image resolution from 100 to 400 dpi. The curves are very similar to each other.

32

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

Table 3 Results of inversions for three slices of Palmottu granite and for different parameterisation of diffusion Slice 2 One diffusion coefficient Two diffusion coefficients

One power factor Two power factors

Slice 5 − 11

2 −1

D = 4.8 × 10 m s R = 1.8 × 10− 2 DPatches = 4.3 × 10− 13 m2 s− 1 DFissures = 10− 9 m2 s− 1 = D0 R = 3 × 10− 3 m = 0.60 R = 1.7 × 10− 2 mPatches = 1.86 mFissures = 0 R = 6.5 × 10− 4

with two variables (two Da or two m, one for the fissures and one for the porous patches), is an improvement but depends on the slice used as reference map. The fitting for slice 2 is very good when using two

Slice 7 − 11

2 −1

D = 3.5 × 10 m s R = 1.9 × 10− 2 DPatches = 9.5 × 10− 14 m2 s− 1 DFissures = 10− 9 m2 s− 1 = D0 R = 4.4 × 10− 3 m = 0.64 R = 1.6 × 10− 2 mPatches = 2.37 mFissures = 0 R = 1.9 × 10− 3

D = 10 × 10− 11 m2 s− 1 R = 2.2 × 10− 2 DPatches = 9 × 10− 12 m2 s− 1 DFissures = 10− 9 m2 s− 1 = D0 R = 1.7 × 10− 2 m = 0.43 R = 2.1 × 10− 2 mPatches = 0.94 mFissures = 0 R = 1.5 × 10− 2

power law factors (R = 6.5 × 10− 4). This proves that according to the pore microstructure, only two parameters are sufficient to correctly fit the experiment. We will return to the problem of representativity

Fig. 11. Comparison between simulated and experimental out-leached curves after inversion. Solid lines represent simulated curves, data points are derived from experiment.

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

of 2D slices used to model an experiment performed in 3D. • The quality of a fitting is influenced by the type of variable. Using Archie's law for the diffusion coefficients always provides a fitting closer to the data than using diffusion coefficients independent of the porosity. This result was observed to be independent of the number of variables and identity of the slice. To help to understand this effect, take the example of a fissure containing a throat. If a unique diffusion coefficient is imposed for all the pixels of a fissure, the constrictivity associated with the throat will not be accounted for by geometry factor G, even if the porosity of the throat is lower than in the rest of the fissure. Imposing a power law relationship between G and Φ (i.e. expression (12)) would decrease the value of G (or increase the constrictivity) near the throat. At pixel scale, our results argue in favour of geometry factors G depending on porosity. Authors that undermine Archie's relation have not used it correctly, i.e. in simulations of diffusion that account for the geometry of a porous network. • A fitting varies according to the 2D slices used as reference maps of a porous network. This last result is globally in agreement with the simulation results obtained in Section 2.5. Although the three slices tested here were obtained from the same sample, they were selected especially because of their differing spatial distributions of porosity (see Section 2.1). These differences were due to the poor representativity of the samples in terms of volume as compared to the mean grain size of the rock. When using the

33

microstructure of slice 7 to fit the experimental data, it is impossible to get a satisfactory result even with two power law exponents. For instance, slice 7 intersects numerous porous patches. The comparison with the simulation obtained from the virtual 3D sample (see Section 2.5) also suggests that the anomalous connectivity of such patches with the border of the slice prevents a correct fitting of the experimental curve. On the other hand, calculations regarding slices 2 and 5 fit accurately the experimental curve when two power law factors are used. These two latter slices (especially slice 2) are much more representative of the porosity pattern observed in the sample, and unlike slice 7, the porous patches are not well connected to the border of the sample. An image of the particle cloud at different times helps understanding how the spatial distribution of porosity may influence the shape of the leaching curve (Fig. 12). The initial portion of a leaching curve (e.g. 0–750 s1/2) corresponds to particles leaching from the fissures and from porous patches directly connected to the border of the sample (Figs. 11 and 12). The remainder of the curve (e.g. 750–3500 s1/2) relates to particles initially located in a reservoir made of porous patches within the sample and connected to the border of the sample by fissures. • The orders of magnitude of the sought diffusion coefficients are realistic. For a single variable, inversion is not precise (Fig. 7a,b). However, the apparent diffusion coefficients are within the classical range of crystalline rocks, i.e. 10− 10–10− 11 m2 s− 1. The values for the most representative slices (slices 2 and 5) yielded from unique Archie coefficients (m = 0.60–

Fig. 12. Particle cloud computed over the network of slice 2 at two times. Black pixels represent one or several particles. Da = 10− 9 m2 s− 1 and Da = 10− 11 m2 s− 1 in fissures and porous patches, respectively. The image is 416 × 416 pixels (including the free fluid around the sample), and the pixel size is 127 × 127 μm2.

34

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

0.64) are very similar to the empirical coefficient found by Parkhomenko (1967) and Sato (1999): m = 0.58. Even if the fitting with a single exponent m is fair when using the geometry of the porous network, the optimal m determined from inversion is comparable to experimental results. Applying Archie's law to the bulk porosity of the sample (Φ = 0.5%, m = 0.58) yields Da = 4.6 × 10− 11 m2 s− 1. This value is equal to the Da found for slices 2 and 5, but is five times less than the value obtained directly from the leaching curve using the Hespe (1971) method: Da = 1.3 × 10− 10 m2 s− 1 (Sardini et al., 2003). For fittings with two variables, fissures and porous patches were analysed separately. Whatever the analysis method, the apparent diffusion coefficients of fissure subset pixels are equal to the free diffusion coefficient (Da =D0, m = 0). This result is logical if we consider a 2D fissure intersecting a 127 × 127 μm2 pixel. The apertures of the fissures in Palmottu granite are micrometric. The overlapping between a fissure and a pixel at this scale is almost a straight stripe, consequently with a tortuosity close to unity. When a unique diffusion coefficient is considered for the porous patches the geometry factors are very low: 1/2300 for slice 2. This is indicative of a very high tortuosity and constrictivity in these porous subsets yielding high values of the m exponent of Archie's law: 1.82 and 2.37 for slices 2 and 5, respectively. Using m = 1.82, the distribution of Da = D0Φ m was recalculated over the porous patches of slice 2 (Fig. 13). The distribution of Da is roughly log normal, with a geometric mean of 8.4 × 10− 12 m2 s− 1 (G = 1/119). We use here the geometric mean by analogy with the best equivalent permeability computed for randomly distributed

permeability fields (Renard and de Marsily, 1997). This value is more reasonable than 1/2300, because it is more consistent with the order of magnitude of Da measured from finely divided materials without sorption (e.g. see Van Loon et al. (2003) for Opalinus Clay from Mont Terri, Switzerland, Melkior et al. (2005) for Bure Argilite, France, and Tits et al. (2003) for cement pastes). Furthermore, a synthesis of diffusivity data within rocks made by Nakashima (1995) shows the existence of two rock types, named type A and B by the author. Type A was related to rock type having pores larger than 25 nm, diffusion coefficients being fitted with m = 0.3, and type B was related to rock type having pores smaller than 25 nm, diffusion coefficient being fitted with m = 2. This “bulk” result is in very good agreement with the present estimation of power-law exponents for microcracks (m = 0) and for altered minerals (m = 1.82). 4. Discussion — conclusions 4.1. Representation of diffusion from 2D to 3D To our knowledge, this study represents the first diffusion experiment in a crystalline rock that has been inverted by conditioning the problem using an accurate spatial representation of the porosity distribution. The Palmottu granite was chosen in order to characterise transport in heterogeneous rock composed of phases with different properties (Siitari-Kauppi et al., 2003). This work highlights the critical role of porous components (altered minerals) with high internal surface areas in very small pores embedded in a conductive

Fig. 13. Distribution of apparent diffusion coefficients in the porous patches of slice 2 from Archie's law with exponent m = 1.86 and a porosity distribution obtained from the 14C-PMMA method.

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

network (fissures). The analytical approach used 2D maps to model the diffusion process, although the experiments were essentially 3D. It is still almost impossible to apply this approach to the same scale in 3D because (i) a representative volume of pore space in crystalline rocks is not yet available, and (ii) specific developments are required to make 3D porosity maps compatible with reasonable computation times. However, in the case of the Palmottu granite, we assumed that the 2D approach was relevant. The 2D representation of the pore network of the Palmottu granite is clearly above the 2D critical percolation threshold, and the 2D pore network presumably forms an “infinite cluster” (Stauffer and Aharony, 1994). An indicator of the internal structure of such a cluster is the fraction of 2D dead-ends, which represents only 5–7% of the total connected porosity network (Fig. 9). This value is very

35

low, and corresponds to well-connected networks. It remains, however, that numerous other crystalline rocks show a 3D 14C-PMMA impregnated pore network not well-connected in 2D. Two typical examples of 3D connected pore networks that are not well-connected in 2D are presented in the left part of Fig. 14. The Le Bourneix granite presents a 3D connected pore network mainly constituted by microporous minerals that are not 2D connected, but that are 3D connected. Although mainly composed of microcracks, the connected pore network of the Charroux tonalite is not 2D connected. The 3D reconstruction of the pore space for these two rocks cannot be dismissed in order to simulate well conditioned matrix diffusion. Conversely, the connected porosities of the Grimsel granodiorite and Palmottu granite (right of Fig. 14) appear to be well-connected in a 2D section.

Fig. 14. Some commonly encountered types of 3D connected pore space of crystalline rocks using the 14C-PMMA method. The pore space of the Charroux tonalite and Le Bourneix granite is not 2D connected, whereas the pore space of the Palmottu granite and Grimsel granodiorite is 2D connected. Scale bars: 1 cm.

36

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37

4.2. Diffusion coefficients and residence time distribution in the matrix The best curve fitting of the experimental outdiffusion data was obtained for the most representative slice (empirical appraisal at glance), with Da = D0 × Φ 1.86 for porous patches and Da = D0 for fissures (D0 = 10 − 9 m 2 s − 1 ). Knowing the porosity distribution in the porous patches, it was then possible to determine the distribution of Da in these porous subsets (Fig. 13). Haggerty et al. (2000, 2001) have pointed out the importance of characterising the distribution of diffusion coefficients in the matrix because they have a strong impact on the interpretation of tracer tests. Our results are encouraging, but they should be used at the scale of in-situ tracer tests. The distribution of matrix diffusion coefficients and the pore network geometry control the residence time distributions. In fact, these distributions are often evoked to explain the tailing effects observed on breakthrough curves, and they have been expressed by analytical functions chosen a priori (Haggerty et al., 2000). Our work proposes an alternative method in upscaling diffusion properties from laboratory to field scale. In order to model transport at the scale of an insitu experiment, the by-product of diffusion coefficient distributions drawn from inversion at the laboratory scale could be compared to the analytical distributions of matrix diffusion coefficients. The obvious benefit of the proposed approach is the insurance that porosity microstructure will be well conditioned on the matrix heterogeneity. Further works will address this issue on rock matrices located near fractures and characterised by non-stationary spatial properties. Acknowledgments We are grateful to Esa Oila for helping us with the experiments and the computer code. The authors wish to thank the two anonymous reviewers for their relevant remarks. Support for this research was provided by a STUK grant (Finnish Authority for Radiation and Nuclear Safety). References Bear, J., Tsang, T.F., de Marsilly, G., 1993. Flow and Contaminant Transport in Fractured Rock. Academic Press, San Diego. Carrera, J., Sanchez-Vila, X., Benet, I., Medina, A., Galarza, G., Guimera, J., 1998. On matrix diffusion: formulations, solution methods and qualitative effects. Hydrogeol. J. 6, 178–190. Coster, M., Chermant, J.L., 1989. Précis d'analyse d'images. Presses du CNRS.

Delay, F., Porel, G., 2003. Inverse modeling in the time domain for solving diffusion in a heterogeneous rock matrix. Geophys. Res. Lett. 30, 1147. Delay, F., Porel, G., Sardini, P., 2002. Modelling diffusion in a heterogeneous matrix with a time domain Lagrangian method and an inversion procedure. C.R. Géosci. 334, 967–973. Haggerty, R., 2002. Matrix diffusion : heavy tailed residence time distributions and their influence on radionuclide retention. Radionuclide Retention in Geological Media, Workshop Proceedings, Oskarshamn, Sweden, pp. 81–90. Haggerty, R., McKenna, S.A., Meigs, L.C., 2000. On the late-time behaviour of tracer test breakthrough curve. Water Resour. Res. 36, 3467–3479. Haggerty, R., Fleming, S.W., Meigs, L.C., McKenna, S.A., 2001. Tracer tests in a fractured dolomite 2. Analysis of mass transfer in single-well injection — withdrawal tests. Water Resour. Res. 37, 1129–1142. Hellmuth, K.H., Siitari-Kauppi, M., Lindberg, A., 1993. Study of porosity and migration pathways in crystalline rocks by impregnation with 14C-polymethylmethacrylate. J. Contam. Hydrol. 13, 403–418. Hespe, E.D., 1971. Leach testing of immobilized radioactive waste solids. In: Turkov, Z.I. (Ed.), Atomic Energy Review, vol. 9, pp. 195–207. Jacob, A., Mazurek, M., Heer, W., 2003. Solute transport in crystalline rocks at Aspö — 2: blind predictions, inverse modelling and lessons learnt from test STT1. J. Contam. Hydrol. 61, 175–190. Johansson, H., Siitari-Kauppi, M., Skalberg, M., Tullborg, E.L., 1998. Diffusion pathways in crystalline rock. Examples from Äspödiorite and fine-grained granite. J. Contam. Hydrol. 35, 41–53. Mazurek, M., Jakob, A., Bossart, P., 2003. Solute transport in crystalline rocks at Aspö — 1: geological basis and model calibration. J. Contam. Hydrol. 61, 157–174. Melkior, T., Yahiaoui, S., Motellier, S., Thoby, D., Tevissen, E., 2005. Cesium sorption and diffusion in Bure mudrock samples. Appl. Clay Sci. 29, 172–186. Müller, G., 1967. Methods in sedimentary petrology. In: Von Engelhardt, W., Füchtbauer, H., Müller, G. (Eds.), Sedimentary Petrology, Part 1, Stuttgart (Schweitzerbart'sche Verlagsbuchhandlung). Nakashima, S., 1995. Diffusivity of ions in pore water as a quantitative basis for rock deformation rate estimates. Tectonophysics 245, 185–203. Neretnieks, I., 1980. Diffusion in the rock matrix: an important factor in radionuclide retardation? J. Geophys. Res. 85, 4379–4397. Oila, E., Sardini, P., Siitari-Kauppi, M., Hellmuth, K.H., 2005. The 14 C-PMMA impregnation method and image analysis as a tool for porosity characterisation of rock-forming minerals. Geological Society London, Special Publication. Petrophysics of Crystalline Rocks, vol. 240, pp. 335–342. Parkhomenko, E.I., 1967. Electrical Properties of Rocks. Plenum Press, New York. Renard, P., de Marsily, G., 1997. Calculating equivalent permeability: a review. Adv. Water Resour. 20, 253–278. Sardini, P., Moreau, E., Sahel, A., Badri, A., Touchard, G., Meunier, A., Beaufort, D., 1997. Two-dimensional fracture and porous network geometry: detection of dead-ends using morphological operations. Acta Stereol. 16, 31–41. Sardini, P., Moreau, E., Sammartino, S., Touchard, G., 1999. Primary mineral connectivity of polyphasic igneous rocks by high-quality digitisation and 2D image analysis. Comput. Geosci. 25, 599–608. Sardini, P., Delay, F., Hellmuth, K.H., Porel, G., Oila, E., 2003. Interpretation of out-diffusion experiments on crystalline rocks using random walk modeling. J. Contam. Hydrol. 61, 339–350. Sato, H., 1999. Matrix diffusion of some simple cations, anions, and neutral species in fractured crystalline rocks. Nucl. Technol. 127, 199–211.

P. Sardini et al. / Journal of Contaminant Hydrology 93 (2007) 21–37 Siitari-Kauppi, M. 2002. Development of 14C-polymethylmethacrylate method for the characterisation of low porosity media. Application to rocks in geological barriers of nuclear waste storage. PhD thesis, Universities of Helsinki and Poitiers. Siitari-Kauppi, M., Marcos, N., Klobes, P., Goebbels, J., Timonen, J., Hellmuth, K.H., 2003. The Palmottu Natural Project. Physical Rock Matrix Characterisation. Geological Survey of Finland, Report YST-118. . Skagius, K., Neretniecks, I., 1986. Porosities and diffusivities of some nonsorbing species in crystalline rocks. Water Resour. Res. 22, 389–398. Stauffer, D., Aharony, A., 1994. Introduction to Percolation Theory. Taylor and Francis, London. Tidwell, V., Meigs, L., Christian-Frear, T., Boney, C., 2000. Effects of spatially heterogeneous porosity on matrix diffusion as investigated by X-ray absorption imaging. J. Contam. Hydrol. 42, 285–302.

37

Tits, J., Jacob, A., Wieland, E., Spieler, P., 2003. Diffusion of tritiated water and 22Na through non-degraded hardened cement pastes. J. Contam. Hydrol. 61, 45–62. Van Loon, L.R., Soler, J.M., Bradbury, M.H., 2003. Diffusion of HTO, 36 Cl and 125I in Opalinus Clay samples from Mont Terri. Effect of confining pressure. J. Contam. Hydrol. 61, 73–83. Vilks, P., Cramer, J.J., Jensen, M., Miller, N.H., Miller, H.G., Stanchell, F.W., 2003. In situ diffusion experiment in granite: phase 1. J. Contam. Hydrol. 61, 191–202. Worthington, P.F., 1993. The uses and abuses of the Archies equations, 1: the formation factor–porosity relationship. J. Appl. Geophys. 30, 215–228.