Inverse modeling of tracer experiments in FEBEX compacted Ca-bentonite

Physics and Chemistry of the Earth 31 (2006) 640–648 www.elsevier.com/locate/pce

Inverse modeling of tracer experiments in FEBEX compacted Ca-bentonite Javier Samper

a,*

, Zhenxue Dai a,1, Jorge Molinero b, M. Garcı´a-Gutie´rrez c, T. Missana c, M. Mingarro c

a

c

Escuela Te´cnica Superior de Ingenieros de Caminos, Canales y Puertos, Campus de Elvin˜a s/n, 15192 La Corun˜a, Spain b Escola Polite´cnica Superior, Universidad de Santiago de Compostela, 27002 Lugo, Spain Dpto. de Impacto Ambiental de la Energı´a, Caracterizacio´n Hidrogeoquı´mica de Emplazamientos, CIEMAT, Edificio 20-A, Avda. Complutense, 22, 28040 Madrid, Spain Received 18 September 2005; received in revised form 28 January 2006 Available online 22 June 2006

Abstract Solute transport parameters of compacted Ca-bentonite used in the FEBEX Project were derived by Garcı´a-Gutie´rrez et al. (2001) from through- and in-diffusion experiments using analytical solutions for their interpretation. Here we expand their work and present the numerical interpretation of diffusion and permeation experiments by solving the inverse transport problem which is formulated as the minimization of a weighted least squares criterion measuring the differences between computed and measured concentration values. The inverse problem is solved with INVERSE-CORE2D, a finite element code which accounts for both dissolved and sorbed concentration data, uses either the Golden section search or Gauss–Newton–Marquardt methods for minimizing the objective function and allows the estimation of transport and retardation parameters such as diffusion coefficient, total and kinematic porosity and distribution coefficients. Diffusion and permeation experiments performed on FEBEX compacted bentonite using tritium, cesium, selenium, and strontium have been effectively interpreted by inverse modeling. Estimated parameters are within the range of reported values for these tracers in bentonites. It has been found that failing to account for the role of sinters may lead to erroneous diffusion coefficients by a factor of 1.4. Possible ways to improve the design of in-diffusion and permeation experiments have been identified. The interpretation of the tritium permeation experiment requires the use of a double-porosity model with mobile porosity of 0.14 for a dry density of 1.18 g/cm3.  2006 Elsevier Ltd. All rights reserved. Keywords: Inverse modeling; Solute transport parameters; FEBEX; Bentonite; Diffusion experiments; Sorption

1. Introduction Predictions of reactive solute transport through compacted bentonites are needed for the performance assessment of the disposal of hazardous wastes such as highlevel nuclear waste (HLW). The FEBEX bentonite was extracted from the Cortijo de Archidona deposit, exploited * Corresponding author. Address: Escuela de Caminos, Campus de Elvin˜a s/n, Universidad de La Corun˜a, 15192 La Corun˜a, Spain. Tel.: +34 981 16 70 00x1433; fax: +34 981 16 71 70. E-mail address: [email protected] (J. Samper). 1 Present address: Earth and Environmental Sciences Division, Los Alamos National Laboratory, New Mexico, USA.

1474-7065/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.pce.2006.04.013

by Minas de Ga´dor, S.A. in Serrata de Nı´jar (Almerı´a, Spain). This bentonite was selected by ENRESA (Empresa Nacional de Residuos Radioactivos, S.A.) prior to the FEBEX project (Huertas et al., 2000) as a suitable material for backfilling and sealing of a HLW repository because it has a very large montmorillonite content, large swelling pressure and sorption capacity, extremely low permeability, acceptable thermal conductivity and ease of compaction for the fabrication of blocks. Bentonite used in the FEBEX project was homogenized to reduce the uncertainties in parameter variability (Huertas et al., 2000). The solution of the inverse problem provides a way whereby measurements of state are used to determine unknown flow and transport parameters by fitting model

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

outputs to measurements. A general framework and description of the coupled inverse problem of water flow and solute transport is given by Mishra and Parker (1989) and Sun (1994). More work on inverse flow and solute and transport in aquifers has been done by Samper et al. (1990), Gailey et al. (1991), Wagner (1992), Xiang et al. (1993), Medina and Carrera (1996) and Poeter and Hill (1997). However, the ability to fully characterize transport parameters of compacted bentonite has not kept pace with the numerical and modeling expertise. Deriving reliable estimates of transport and sorption parameters of the FEBEX bentonite is one of the objectives of the FEBEX Project. Experimental characterization of diffusive transport of radionuclides in FEBEX bentonite was carried out by means of through- and in-diffusion experiments (Garcı´a-Gutie´rrez et al., 2001). Application of the inverse approach to estimate various flow and transport parameters of the FEBEX bentonite are presented by Dai and Samper (1999). Solute transport parameters of compacted Ca-bentonite used in the FEBEX Project were presented by Garcı´a-Gutie´rrez et al. (2001). These parameters were derived from through- and in-diffusion experiments which were interpreted by means of analytical solutions. Here we expand on the latter approach and provide results not only for diffusion but also for flow-through or permeation experiments. Numerical interpretation of laboratory experiments is achieved by solving the inverse problem of solute transport. The inverse approach overcomes the limitations of analytical solutions by accounting for: (1) Measurement errors; (2) The role of the sinters; (3) Time-varying boundary conditions in permeation experiments; and (4) Transient concentrations in the cell and total concentration data at the end of in-diffusion experiments. The inverse analysis has been applied to the numerical interpretation of through-diffusion, in-diffusion and permeation experiments performed on FEBEX Ca-bentonite using tritium, strontium, cesium and selenium. The following solute transport parameters are estimated: molecular diffusion coefficient, dispersivity, distribution coefficient, Kd, and kinematic as well total porosities. 2. Direct problem of solute transport and sorption The equation for the transient transport of a sorbing solute through a saturated porous medium under steady-state flow is given by (Bear, 1972): r  ð/DrCÞ  qrC þ rðC   CÞ ¼ /R

oC ot

ð1Þ

where / is porosity, C is solute concentration; q is volumetric water flux; r is a fluid sink/source term having a concentration C*, t is time, and R is the retardation coefficient defined as R¼1þ

qK d /

ð2Þ

641

where Kd is the distribution coefficient and q the bulk density. Tracer concentrations are extremely small. The largest concentrations are 2.69 · 109 M for cesium, 5.4 · 109 M for selenium and 109 M for strontium. For such low concentrations all sorption isotherms are linear. The dispersion tensor, D, in (1) is given by D ¼ De I þ Dh

ð3Þ

and includes the hydrodynamic or mechanical dispersion tensor, Dh, and the molecular diffusion term. Since clay samples were prepared by compacting powder FEBEX bentonite which had been previously homogenized, molecular diffusion is expected to be isotropic. Therefore, the diffusion term can be written as DeI, where De is the effective molecular diffusion coefficient and I is the identity tensor. In one-dimensional solute transport the mechanical dispersion, Dh, is equal to aq where a is dispersivity. The effective diffusion coefficient in porous media, De, is usually expressed as: De ¼ /Dp

ð4Þ

where Dp is the molecular diffusion coefficient in the pore space which in turn is given by Dp ¼ D0 s

ð5Þ

where D0 is the molecular diffusion coefficient in pure water and s is the tortuosity of the medium which is assumed here to be equal to /1/3 where / is porosity. The apparent diffusion coefficient, Da, takes into account the combined effect of diffusion and sorption and is defined as Da ¼

De / þ qK d

ð6Þ

Solution of Eq. (1) requires appropriate initial conditions: Cjt¼0 ¼ C 0

ð7Þ

and boundary conditions: /DrC  nC1 ¼ bðC 1  CÞ þ F 0D

ð8Þ

where n is a unit vector normal to the boundary C1 pointing outwards, C0 is initial concentration, C1 is a prescribed concentration; b is a parameter controlling the type of boundary condition. For a Neumann condition, b = 0 and solute flux through the boundary is equal to the dispersive flux, F 0D . For a Dirichlet condition, b = 1, F 0D ¼0 and the concentration is prescribed to be equal to C1. At water inflow boundaries b = q Æ n and F 0D ¼0. 3. Inverse problem of solute transport and sorption The essence of the inverse problem lies on deriving optimum parameter estimates from known concentration data. Optimum parameters are those which minimize an objective function measuring the difference between measured and computed concentrations. Our formulation of the inverse problem is based on generalized least squares criterion (Sun, 1994; Dai and

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

Samper, 2004). Let P = (p1,p2,p3, . . . ,pM) be the vector of M unknown parameters. The least squares criterion, E(p), can be expressed as, EðpÞ ¼

NE X

W i Ei ðpÞ

ð9Þ

i¼1

where i ¼ 1; NE denotes different types of data, i = 1 for dissolved concentration; i = 2 for total concentration; i = 3 for prior information on model parameters, Wi is the weighting coefficient of the ith generalized least-squares criterion Ei(p) which is defined as: Li X w2li r2li ðpÞ Ei ðpÞ ¼ l¼1

rli ðpÞ ¼ uil ðpÞ  F il where uil ðpÞ is the computed value of the ith variable at the lth observation point; F il are measured values; Li is the number of observations, either in space or in time, for the ith type of data; and rli is the residual or difference between model outcome and measurement; wli is the weighting coefficient for the lth measurement of the ith type of data. Its values depend on the accuracy of observations. If some data are judged to be unreliable, they are assigned small weights in order to prevent their pernicious effect on the optimization process. Eq. (9) is a weighted multiobjective optimization criterion. Weighting coefficients Wi are computed from (Carrera and Neuman, 1986; Dai and Samper, 2004): Wi ¼

W 0i Ei ðpÞ Li

i ¼ 1; NE

ð10Þ

where W0i are user-defined dimensionless initial weighting coefficients for different types of observation data. While coefficients W0i are dimensionless, Wi have dimensions which are reciprocal of those of Ei(p). Generally Ei(p) has units equal to the square of the unit used for each type of data (i.e., (mol/l)2 for dissolved concentrations). Weights Wi are updated automatically during the iterative optimization process according to Eq. (10) whenever NE > 1. Weights wli for different observation points in (6) must be specified in advance. In a statistical framework, these coefficients are inverse standard deviations expressing the uncertainty of measured data. Such standard deviations can be derived from measurement or analytical errors, rm, or in the case of time series from random fluctuations about a trend. For a given type of data (i.e., for a fixed i), wli for all l = 1, 2, . . . , Li are all equal to r1 m whenever all data values are equally reliable. Minimization of the expression in Eq. (9) can achieved by different methods such as: golden section search (Sun, 1994), Newton, Gauss–Newton, Gauss–Newton–Levenberg–Marquardt (Yeh, 1986; Sun, 1994), conjugate gradient (Carrera and Neuman, 1986), quasi-Newton, and simulated annealing. In this paper, golden section search and Gauss–Newton–Levenberg–Marquardt methods have

been used to solve the inverse problem. Combining these two inverse methods with the forward modeling code CORE2D (Samper et al., 2000; Xu et al., 1999), an inverse code INVERSE-CORE2D (Dai and Samper, 1999) has been developed. The code can estimate any flow and transport parameter in the code CORE2D and provides the statistical measures of goodness-of-fit as well as parameter uncertainty by computing the covariance and correlation matrices, eigenvalues and approximate confidence intervals (see Dai and Samper, 2004). 4. Description of transport experiments Several types of diffusion experiments (Fig. 1) were carried out by Garcı´a-Gutie´rrez et al. (2001) to estimate transport and sorption parameters of FEBEX bentonite, including total and kinematic porosity, diffusion coefficients (effective and apparent) and distribution coefficient. These experiments were interpreted by means of available analytical solutions by Garcı´a-Gutie´rrez et al. (2001). The conditions of some experiments are not taken into account by such solutions. The use of numerical models for experiment interpretation overcomes these limitations of analytical solutions. 4.1. Through-diffusion experiments Through-diffusion experiments were performed using tritium (experiments TD-HTO1, TD-HTO2 and TDHTO3) and strontium (experiments TD-Sr5, TD-Sr6, TD-Sr7 and TD-Sr8). Experiments were performed on bentonite plugs of 50 mm of diameter, with lengths ranging from 5.3 mm for

PERMEATION

Pressure

642

TRACER

CLAY

C

THROUGH DIFFUSION

TRACER IN

CLAY

TRACER OUT

CLAY

TRACER IN

IN DIFFUSION

TRACER IN

Fig. 1. Methods used for determining sorption and diffusion parameters (from Garcı´a-Gutie´rrez et al., 2001).

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

tritium to 8.3 mm for strontium (Fig. 1). They were immersed in water for 4 weeks achieving a final density of 1.18 g/cm3. The clay plug was located in contact with two cells. To ensure a uniform diffusion through the whole cross-section of the plug, porous stainless-steel sinters manufactured by Mott were located at both sides of the plug. According to technical specifications sinters are 3.35 mm long and have a porosity of 0.5. A known mass of tracer was added initially to one of the cells. The evolution of tracer activity was measured in both cells at regular intervals. Experiments lasted around 54 days for tritium and nearly 200 days for strontium. The model domain of through-diffusion experiments includes the following subdomains: the upstream and downstream cells, two sinters and the clay plug. Finite elements representing the cells were assigned a porosity of 1 and a large-enough diffusion coefficient in order to ensure that all nodes within the cell have the same activity. For the sinters the effective diffusion was computed using Eqs. (4) and (5) with a porosity of 0.5. The initial activity in the upstream cell was taken equal to the initial prescribed tracer activity. Initial activity is equal to zero in the rest of the subdomains. Two finite element grids were used for the numerical simulation of through diffusion experiments. That of tritium (TD-HTO1, TD-HTO2 and TD-HTO3) includes 200 triangular elements and 202 nodes. The grid for strontium experiments (TD-Sr5, TD-Sr6, TD-Sr7 and TD-Sr8) has 114 triangular elements and 116 nodes. The simulation time horizon covers the expected duration of the experiments. For TD-HTO1, TD-HTO2 and TD-HTO3 the total time is 54 days, which was discretized with 800 time increments. For strontium experiments the total experiment time of 193 days was divided into 1800 time increments. Before running the inverse model, some forward runs were made to test and ensure the numerical stability of the solutions.

643

Two kinds of numerical models were used in these experiments, one for cesium (In-Cs-9 and In-Cs-10) with 600 triangular elements and 602 nodes, and the other for selenium (In-Se-17, In-Se-18) with 200 triangular elements and 202 nodes. The simulation time for both tracers (cesium and selenium) is 439 days and it is discretized with 630 time increments. 4.3. Permeation experiments A high-pressure stainless steel cell was used for permeation experiments. The clay plug, 5.3 mm in length and 50 mm of diameter was sealed between two 3.35 mm thick stainless steel sinters (Fig. 1). The granite synthetic groundwater was injected into the sample with a low driving pressure (6 bar) ensured by a HPLC pump with pressure and flow regulation (the flow rate in the cells ranged from 2.7 to 4.0 ml/day). The final density is 1.18 g/cm3. All the experiments were performed under oxic conditions. After 4 weeks of sample saturation, a small pulse of tritium tracer (500 ll or 183,000 counts) was injected into the fluid pathway. The breakthrough curve of the tracer at the outlet was monitored continuously by a fraction collector. Single and double porosity models were used to model permeation experiments. The former is simulated with a 1D grid while the latter requires a 2-D finite element grid containing a mobile part where water can flow and an immobile part where only molecular diffusion can take place. Since the tracer was injected as a non-instantaneous pulse, the duration and distribution of the pulse is uncertain. Only the flow rate and the total tracer mass injected are known. Tracer injection was modeled by using a time-dependent pulse function. A trial and error method was used to estimate such pulse function (Fig. 2) by keeping constant the total mass of tracer injected. Initial activities in all of the nodes are zero.

4.2. In-diffusion experiments 1000

Permeation 1 (P-1d) (Double-porosity) Time function 1 Obj.=185.6 Time function 2 Obj.=72.0 Time function 3 Obj.=12.2

800

Activity (Mcount/l)

In-diffusion (ID) experiments consist on the immersion of a clay sample of 38 mm of diameter and 25 or 50 mm of length into a cell containing granite synthetic groundwater until saturation. Then a tracer, which is initially not present either on the clay or on the sinters, is added to the solution and allowed to diffuse into clay samples from both sides. The gradient in tracer activities induces a diffusive flux from the cell into clay samples (Fig. 1). The decrease in tracer activity in the immersion cell was monitored continuously. At the end of the experiment clay plugs were sliced to obtain the activity profile along clay samples. In-diffusion experiments were carried using cesium and selenium. Se(IV) was added as a sodium selenite SeO3Na2. Tests were performed under oxic conditions. The model domain for in-diffusion experiments contains five parts: two cells, two sinters and a clay plug. Initial activities in the two cells are equal to the initial tracer activity while in the sinters and the clay plug they are equal to zero.

600

400

200

0 0.0

0.5

1.0

1.5

2.0

Time (days) Fig. 2. Time functions used to define the tracer pulse of the inflow water. Also shown are the values of the objective function for each time function. One can see the large sensitivity of the fit to the selected time function.

644

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

per, 1999). Results of the numerical interpretation of TD, IN and P experiments are shown in Table 1. This table contains also the parameters obtained by Garcı´a-Gutie´rrez et al. (2001) by using analytical methods. It should be noticed that not all experiments could be interpreted analytically. Fig. 3. Finite element mesh for the double-porosity model.

The 1-D single porosity model (P-1) has 200 triangular elements and 202 nodes. This model proved to be unable to match observed data and then a 2-D double porosity model (P-1d) was used which was discretized with 560 triangular elements and 324 nodes. Fig. 3 illustrates the finite element mesh used for the double-porosity model. The total length of the system is 1.2 cm, including two sinters at both sides (2 · 0.335 cm long) and the bentonite sample in the middle (0.53 cm long). The width is 0.2 cm. Based on a prior estimate of mobile porosity of about 20%, the width of the mobile part (bottom row) is 0.04 cm while that of the immobile part (the upper seven rows) is 0.16 cm. The simulation time is 3.78 days. It is discretized with 88 time periods with gradually increasing time increments. 5. Interpretation by inverse modeling Diffusion experiments were interpreted numerically using the inverse code INVERSE-CORE2D (Dai and Sam-

5.1. Through-diffusion experiments Several sets of through-diffusion (TD) experiments were performed on clay samples at dry densities ranging from 1.18 to 1.65 g/cm3 using tritiated water (tritium) and strontium. Tritium was used in three TD experiments which lasted about 50 days, a time long enough for the experiment to attain steady activities along the samples. Tritium experiments were used to test the estimation algorithms and to explore the role of the sinters. Our results indicate that the sinters should be taken into account for a proper interpretation of the experiments. The results in Table 1 clearly demonstrate that bentonite diffusion coefficients obtained from analytical solutions underestimate the true diffusion coefficient by a factor of 1.4. To demonstrate the accuracy of the numerical results, experimental data were interpreted numerically with a model without sinters. In this case, inverse estimates of diffusion coefficients coincide with those obtained with analytical methods. Therefore, failing to account for the sinters can lead to the underestimation of diffusion coeffi-

Table 1 Summary of interpretation of CIEMAT diffusion experiments performed on samples of compacted FEBEX bentonitea Tracer

Experiments and bentonite density (g/ cm3)

Numerical /

De (m /s)

TD-HTO1

1.18

TD-HTO2

1.18

TD-HTO3

1.18

P-1

1.18

0.56b 0.5 0.56b 0.45 0.56b 0.44 0.52/0.14

1.41 · 1010 1.41 · 1010 1.24 · 1010 1.24 · 1010 1.18 · 1010 1.17 · 1010 3.5 · 109

Strontium

TD-Sr5 TD-Sr6 TD-Sr7 TD-Sr8

1.39 1.39 1.65 1.65

0.48b 0.48b 0.38b 0.38b

1.10 · 109 1.81 · 109 1.25 · 109 9.70 · 1010

6.74 · 1013 7.13 · 1013 5.47 · 1013 5.28 · 1013

Cesium

ID-Cs-9

1.65

ID-Cs-10A

1.57

ID-Cs-10B

1.57

0.42b 0.49 0.42b 0.49 0.42b 0.45

8.03 · 1010 8.05 · 1010 3.3 · 1010 3.44 · 1010 7.57 · 1010 8.81 · 1010

3.34 · 1013 3.32 · 1013 1.37 · 1013 1.43 · 1013 3.21 · 1013 3.39 · 1013

ID-Se-17

1.57

ID-Se-18

1.57

0.42b 0.21 0.42b 0.25

4.35 · 1013 1.73 · 1013 1.39 · 1013 1.23 · 1013

6.51 · 1014 3.4 · 1014 2.29 · 1014 2.23 · 1014

Tritium

Selenium

a

Analytical 2

2

Da (m /s)

Kd (ml/g)

De (m2/s)

Da (m2/s)

Kd (ml/g)

8.81 · 1012 3.44 · 1012 3.99 · 1012 5.61 · 1012

1410 1200 1550 1300

3.07 · 1013

823 ± 121

1.03 · 1010 8.89 · 1011 9.80 · 1011

1176 1823 1383 1114 925.1 931.8 924.1 924.4 907.9 999.6 2.41 1.87 2.17 2.02

2.51 · 1013 3.51 · 1013 6.25 · 1014

3.35 ± 1.5

Analytical results were obtained by Garcı´a-Gutie´rrez et al. (2001). Here / denotes diffusion-accessible porosity; De is effective diffusion coefficient; Da is apparent diffusion coefficient and Kd is distribution coefficient in ml/g. b Parameter is fixed and therefore not estimated.

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

cients. Based on this finding, the sinters were considered for the interpretation of all the diffusion experiments. In general, the numerical solution matches observed tracer data in both cells (Fig. 4). Effective diffusion coefficients for tritium range from 1.17 to 1.41 · 1010 m2/s which are within the range of published values for this tracer in bentonites (Yu and Neretnieks, 1996). The 95% confidence interval of the estimated effective diffusion coefficient of TD-HTO1 is rather small: (1.34 · 1010, 1.48 · 1010) m2/ s. The best fit to TD tritium data is obtained with diffusive porosities ranging from 0.44 to 0.5 (see Table 1) which are slightly smaller than the total porosity. Tortuosity factor ranges from 0.76 to 0.79. Estimates of porosity and diffusion (D0), however, are strongly correlated as attested by their large computed correlation coefficients which are greater than 0.95. Given the large uncertainties in porosity and diffusion coefficient, the value of porosity was fixed. Several runs were performed using porosities ranging from 0.45 to 0.56. These runs led to similar fits but different estimates of D0. Estimates of

4.0

Through-Diffusion 1 (TD-HTO1) Observations (IN) Numerical results (IN) Observations (OUT) Numerical results (OUT)

3.5

Activity (Mcount/L)

3.0 2.5 2.0 1.5 1.0 0.5

645

porosity and D0 always provide a similar values of the effective diffusion coefficient. This means that tritium TD experiments only allow the unambiguous estimation of the effective diffusion coefficient. Four TD experiments were performed using strontium: TD-Sr5 and TD-Sr6 with a dry density of 1.18 and TDSr7 and TD-Sr8 with a density of 1.39 g/cm3. This tracer has a large distribution coefficient. In fact, the tracer does not reach the downstream cell after almost 200 days while it has been fully flushed from the upstream cell after that time. These experiments were interpreted in two groups. For TD-Sr5 and TD-Sr6 porosity was fixed to a value of 0.48 that corresponds to a dry density of 1.39 g/cm3. For the second group (TD-Sr7 and TD-Sr8) porosity was fixed to 0.38 that corresponds to a dry density of 1.65 g/cm3. Kd and D0 were estimated in all four experiments. The best fit is shown in Fig. 4. Effective diffusion coefficients De range from 9.7 · 1010 to 1.81 · 109 m2/s with a mean value of 1.28 · 109 m2/s. Estimated values of Kd show also a wide range (from 1114 to 1823 cm3/g) with a mean of 1374 cm3/g. Estimated results of the apparent diffusion Da show much less scatter (from 5.28 · 1013 to 7.13 · 1013 m2/s) with a mean value of 6.15 · 1013 m2/s. Estimation errors and confidence intervals are similar for all four experiments. For experiment TD-Sr5, the 95% confidence intervals for the estimates of effective diffusion and Kd are (5.47 · 1010, 1.65 · 109) m2/s and (786.6, 1566) cm3/g, respectively. These values are within the range of published values for this tracer in bentonites. It should be noticed that in this case the estimates of effective diffusion and Kd show consistently a significant negative correlation of about 0.5. The inverse analysis indicates that Kd is the parameter having the largest uncertainty (estimation error).

0.0 0

10

20

30

40

50

60

Time (days)

25

Activity (Mcount/L)

Several sets of ID experiments were performed on clay samples having a total porosity of 0.42 which amounts to a dry density of 1.57 g/cm3 using radioactive cesium (137Cs) and selenium SeO2þ 3 . Tracer activity in the cell was monitored for full duration of the experiment. At the end of the experiment, clay samples were sliced and total tracer activity was measured in each slice. This activity includes the activity in the liquid and solid phases, in Mcount/per gram of sample (clay + water). In order to fit the transient liquid activity in the cell and the final total activity (including the liquid and solid phases) in the clay at the end of experiments, it is necessary to translate modeling results from liquid activity to total activity by means of

Through-Diffusion 5 (TD-Sr5) Observations (IN) Numerical results (IN) Observations (OUT) Numerical results (OUT)

20

15

10

5

0 0

50

100

150

5.2. In-diffusion experiments

200

Time (days) Fig. 4. Numerical interpretation of through diffusion experiments: tritium TD-HTO1 (top) and strontium TD-Sr5 (bottom). Upstream (IN cell) and downstream (OUT cell) cell activities are shown.

CT ¼ CL

/ þ qd K d /qw þ qd

ð11Þ

where CT = total activity (Mcount/g), CL = liquid activity (Mcount/cm3), qd = bulk dry density of the porous medium (g/cm3), and qw = density of water (g/cm3). Therefore,

646

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

inverse modeling of ID experiments uses two types of data: transient liquid activity in the cell and final total activity in the clay. Two cesium experiments were available (ID-Cs-9 and ID-Cs-10). In one of them, the behavior of the sample is not symmetric, i.e., the tracer diffuses on one side faster than on the other. In this case, data from each side of the sample were interpreted separately (tests ID-Cs-10A and ID-Cs-10B). The experiments were interpreted in two stages. First, the porosity was fixed to a value of 0.42 while Kd and D0 were estimated. In the second stage, porosity was also estimated, resulting in an excellent fit to measured final concentration data. Some discrepancies are observed in the time evolution of tracer activity in the cell which are attributed to uncertainties in the initial activity and the lack of well-mixed conditions in the cell. A detailed sensitivity analysis was performed to evaluate the effects of several sources of uncertainty including: (1) the initial concentration in the sinters, (2) the appropriate boundary condition for the cell, (3) the initial concentration in the cell, and (4) the effective volume of the cell. The results of this analysis indicate that estimated parameters are very sensitive to the properties and conditions prevailing at the cell. Therefore, for the proper interpretation of this type of experiment we suggest that: (1) only one clay sample should be located in each cell, (2) the initial concentration should be measured as accurately as possible, and (3) spe-

cial care should be taken in measuring the activities of the slices located near the faces of the clay samples. The best fit is obtained with diffusive porosities slightly larger (0.49) than total porosity (0.42). This result could be caused by Cs diffusion through clay interlayers. Diffusive porosities range from 0.45 to 0.49. Effective diffusion coefficients De range from 3.44 · 1010 to 8.81 · 1010 m2/s. Estimated values of Kd range from 924.4 to 999.6 cm3/g. Apparent diffusion coefficients range from 1.43 · 1013 to 3.39 · 1013 m2/s. These values are within the range of published values for this tracer in bentonites. Numerical solutions fit perfectly measured data (Fig. 5). Two selenium (in the form of SeO2þ 3 ) experiments were performed on clay samples of 2.5 cm, half the length of the cesium cells. Again, the experiments were interpreted in two stages. First, the porosity was fixed to a value of 0.42 while Kd and D0 were estimated. In the second stage, porosity was also estimated, resulting in an excellent fit to both final concentration data and time evolution of tracer activity in the cell. The best fit is obtained with diffusive porosities smaller than total porosity (Fig. 5). Diffusive porosities are equal to 0.21 and 0.25. Effective diffusion coefficients De range from 1.23 · 1010 to 1.72 · 1010 m2/s. Estimated values of Kd range from 1.87 to 2.02 cm3/g. Apparent diffusion coefficients range from 2.23 · 1013 to 3.40 · 1013 m2/s. Garcı´a-Gutie´rrez et al. (2001) reported difficulties in obtaining a unique fit to sele-

15

In-Diffusion-Cs-9 Measurements Model fit

0.20

Activity (Mcount/l)

Activity (Mcount/g)

0.25

0.15 0.10 0.05

In-Diffusion-Cs-9 Measurements Model fit

10

5

0 0.00 0.0

0.1

0.2

0.3

0.4

0.5

0

100

Distance (dm)

200

300

400

500

Time (days) 30

In-Diffusion-Se-18 Measurements Model fit

0.003

25

Activity (Mcount/g)

Activity (Mcount/g)

0.004

0.002

0.001

0.000

In-Diffusion-Se-18 Measurements Model fit

20 15 10 5 0

0.00

0.05

0.10

0.15

0.20

Distance (dm)

0.25

0.30

0

100

200

300

400

500

Time (days)

Fig. 5. Numerical interpretation of in diffusion cesium ID-Cs-9 (top) and selenium ID-Se-18 (bottom) experiments. (Left figures show the final total activities in the clay sample and figures on the right show the concentration evolution in the cell.)

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

nium IN diffusion experiments that were attributed to chemical reactions. Selenium could be present in several oxidation states and form different complexes with different properties. Although Se(IV) oxidation into Se(VI) seems unlikely, Se(IV) could be reduced by Fe (II) or organic matter inside the clay sample. 5.3. Permeation experiments Contrary to in- and through-diffusion experiments, permeation experiments involve advective and dispersive transport. These experiments are intended to provide information on kinematic porosity which measures the volume of well-interconnected water-flowing pores. Permeation experiments were performed using tritiated water. Following the parsimony principle, the classic single-porosity model was tested first. For the permeation experiment a convergent but suboptimal solution is obtained with a porosity of 0.35 (Fig. 6). Given the limitations of the single-porosity model, a more complex double-porosity model was considered. In our double-porosity model the domain is divided into two domains. Let fm and fim be the fractions of the total

25

PERMEATION 1 (P-1) Measurements Model fit

Activity (Mcount/l)

20

15

10

647

volume of the sample occupied by the mobile and immobile domains, respectively. From prior information we set fm = 0.2 and fim = 0.8. In our formulation of double porosity, model porosities for mobile and immobile zones are not global porosities, but relative porosities. Let /m and /im be the relative porosities of the mobile and immobile zones which are defined as the ratio of void volume to total volume of each domain. The total porosity of the sample is given by / = /mfm + /im fim where /mfm is the overall porosity of the mobile domain and /imfim represents the overall porosity of the immobile domain. According to our experience, this way of parameterizing double porosity is not affected by the geometry adopted in the model. Therefore, the values /mfm and /imfim are independent of the model geometry. The inverse double-porosity model leads to an excellent fit of the activity breakthrough curves (Fig. 6) with relative porosities /m = 0.73 and /im = 0.65 for which the overall porosity of the mobile domain is 0.14 and that of the immobile domain is 0.52. It should be noticed that the estimated total porosity (0.66) is slightly greater than total volumetric porosity (0.56). This could be caused by tritium accumulation in clay minerals which according to Kalinichenko et al. (2002) takes place under room conditions by a two-stage isotope exchange. Tritium ions first migrate from the solution to the bound layers and then substitute the protons of the structural hydroxyls. The effective diffusion coefficient in the mobile phase is 3.5 · 109 m2/s while that of the immobile phase is almost three orders of magnitude smaller. The longitudinal dispersivity for the mobile phase was set equal to one-tenth of the length of the clay sample and was not estimated. 5.4. Comparison of numerical estimations with analytical solutions

5

0 0

1

2

3

4

Time (days)

25

Permeation 1 (P-1d) (Double-porosity) Measurements Model fit

Activity (Mcount/l)

20

15

10

5

0 0

1

2

3

4

Time (days)

Fig. 6. Numerical interpretation of tritium permeation experiment. Computed activities (lines) and measured data (symbols) are shown for single-porosity (top) and double-porosity (bottom) models.

Standard analytical methods (see Garcı´a-Gutie´rrez et al., 2001) for the interpretation of TD and ID experiments do not take into account the sinters. Such a limitation of analytical methods can lead to erroneous estimates of diffusion coefficients. Inverse modeling of tritium TD experiments has revealed that for the experimental conditions of FEBEX bentonite failing to account for the sinters can lead to underestimation of effective diffusion coefficients by a factor of 1.4. Analytical methods for the interpretation of ID experiments rely entirely on activity data measured at the end of the experiment. They make no use of data measured at the immersion cell and allow only the estimation of apparent diffusion coefficient and /R (see Eq. (2)). Inverse modeling, on the contrary, allows the estimation of transport and sorption parameters by using simultaneously both final total activities and time evolution of activities in the immersion cells. Numerical methods are better suited than analytical methods especially for the interpretation of permeation experiments which may require time-varying pulse functions and double porosity models for bentonite.

648

J. Samper et al. / Physics and Chemistry of the Earth 31 (2006) 640–648

Limitations of analytical methods may lead to biased estimates of diffusion coefficients. According to the comparison shown in Table 1, such a bias is about a factor of 1.4 for the effective diffusion coefficient of tritium and an order of magnitude for the apparent diffusion coefficient of strontium. However, the biasness is rather small for the apparent diffusion and the distribution coefficient of Cs and Se (Table 1). 6. Conclusions Diffusion and permeation experiments performed on FEBEX compacted bentonite using tritium, cesium, selenium, and strontium have been effectively interpreted by inverse modeling. Contrary to standard analytical methods, numerical inverse modeling allows an efficient, fast and optimum interpretation of these experiments. It has been found that failing to account for some experimental conditions may lead to erroneous diffusion coefficients which for the experimental conditions of FEBEX bentonite could deviate from true values by a factor of 1.4. Inverse modeling of permeation experiments allows the use of time-varying pulse functions and double porosity models for bentonite. Numerical methods have proved useful for indicating identifiability problems and suggesting possible improvements in experimental conditions such as measuring directly the effective diffusion coefficient of porous sinters. Possible improvements for in-diffusion experiments include: (1) Determining accurately the initial tracer concentration and the effective volume of the cells; and (2) Obtaining detailed and accurate concentration profiles, especially near the faces of the clay samples, and using only one clay sample in each cell. In our work the porosity and the effective diffusion of the sinters were taken from manufacturer technical specifications. However, it is recommended to determine these parameters experimentally by means of diffusion experiments to reduce unneeded uncertainties. The duration of the tracer pulse should be much smaller than the duration of the permeation experiment. Otherwise, the concentration of the inflow water should be carefully measured. Similar to other clays, FEBEX bentonite exhibits a complex porosity structure. Two types of porosities can be distinguished: a mobile porosity along which porewater flows and an immobile porosity where water does not flow while only molecular diffusion takes place. The results of the tritium permeation experiment performed on FEBEX bentonite indicate that mobile porosity is 0.14 for a dry density of 1.18 g/cm3. Acknowledgement This work was supported by the Spanish Nuclear Waste Company (ENRESA) within the framework of the FEBEX Research Project through Research Grants signed with the University of La Corun˜a (Contracts 70323 and 770045).

The FEBEX Project as a whole is funded by the EC (Projects FI4W-CT95-0008 and FIKW-CT-2000-0016) within the Nuclear Fission Safety Programme. We gratefully acknowledge the assistance of John Apps of LBNL. We thank the two anonymous reviewers for their thoughtful reviews of the paper which contributed to improve it. References Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York, 764 pp. Carrera, J., Neuman, S.P., 1986. Estimation of aquifer parameters under steady-state and transient conditions: 1. Background and statistical framework. Water Resour. Res. 22 (2), 199–210. Dai, Z., Samper, J., 1999. INVERSE-CORE2D: A code for the inverse problem of water flow and reactive solute transport. Users Manual, Version 0, University of La Corun˜a, 240 pp. Dai, Z., Samper, J., 2004. Inverse problem of multicomponent reactive chemical transport in porous media. Formulat. Appl. Water Resour. Res. 40, W07407. doi:10.1029/2004WR003248. Gailey, R.M., Gorelick, S.M., Crowe, A.S., 1991. Coupled process parameter estimation and prediction uncertainty using hydraulic head and concentration data. Adv. Water Resour. 14 (5), 301–314. Garcı´a-Gutie´rrez, M., Missana, T., Mingarro, M., Samper, J., Dai, Z., Molinero, J., 2001. Solute transport properties of compacted Ca-bentonite used in FEBEX Project. J. Contam. Hydrol. 47 (3), 127–137. Huertas, F., Fuentes-Cantillana, J.L., Jullien, F., Rivas, P., Linares, J., Farin˜a, P., Ghoreychi, M., Jockwer, N., Kickmaier, W., Martı´nez, M.A., Samper, J., Alonso E., Elorza, F.J., 2000. Full-scale engineered barriers experiment for a deep geological repository for high-level radioactive waste in crystalline host rock (FEBEX project). Final report, EUR 19147, European Communities, 362 pp. Kalinichenko, E.A., Pushkarova, R.A., Fenoll Hach-Alı´, P., Lo´pezGalindo, A., 2002. Tritium accumulation in structures of clay minerals. Clay Miner. 37 (1), 497–508. Medina, A., Carrera, J., 1996. Coupled estimation of flow and solute transport parameters. Water Resour. Res. 32 (10), 3063–3076. Mishra, S., Parker, S.C., 1989. Parameter estimation for coupled unsaturated flow and transport. Water Resour. Res. 25 (3), 385–396. Poeter, E.P., Hill, M.C., 1997. Inverse models: a necessary next step in groundwater modeling. Ground Water 35 (2), 250–260. Samper, F.J., Carrera, J., Galarza, G., Medina, A., 1990. Application of an automatic calibration technique to modeling an alluvial aquifer, ModelCARE 90: Calibration and reliability in groundwater modeling, IAHS Publ. no. 195, pp. 87–95. Samper, F.J., Juncosa, R., Delgado, J., Montenegro, L., 2000. CORE2D: A code for non-isothermal water flow and reactive solute transport, Users Manual version 2, Universidad de La Corun˜a, Techn. Publ. ENRESA, 6/2000, 131 pp. Sun, N.Z., 1994. Inverse Problems in Groundwater Modeling. Kluwer Academic Publishers., Netherlands, 364 pp. Wagner, B.J., 1992. Simultaneous parameter estimation and contamination source characterization for coupled groundwater flow and contaminant transport modeling. J. Hydrol. 135, 275–303. Xiang, Y., Sykes, J.F., Thompson, N.R., 1993. Composite L1 parameter estimator for model fitting in groundwater flow and solute transport simulations. Water Resour. Res. 29 (6), 1661–1673. Xu, T., Samper, J., Ayora, C., Manzano, M., Custodio, E., 1999. Modeling of non-isothermal multicomponent reactive transport in field scale porous media flow systems. J. Hydrol. 214, 144–164. Yeh, W.W.-G., 1986. Review of parameter identification procedures in ground-water hydrology: the inverse problem. Water Resour. Res. 22 (2), 95–108. Yu, J.W., Neretnieks, I., 1996. Diffusion and sorption properties of radionuclides in compacted bentonite, SKB, TR-97-12.