The need for underground repositories has given rise to a number of studies of flow and transport in low-permeability rocks. To isolate the contaminants from the environment, crystalline rocks or clay are considered the most suitable. These media often have large-scale fractures and the safety calculations for the so-called ‘far field’ have to take into account the up-scaling of the flow properties [10] and transport in the fracture networks [1,9]. At the scale of the ‘near field’, it is assumed that the main mechanism affecting potential leaks is diffusion through the rock matrix, which is always heterogeneous. In crystalline rocks, the matrix consists of different grains more or less weathered according to their mineralogy. Micro-cracks and open inter-grain joints may also exist. This results in heterogeneous pore-space patterns and therefore variable pore diffusion coefficients. In clays, the damaged zone around the wells and galleries is large enough to contain unevenly distributed areas of open and connected porosity. Diffusion velocities can be measured either in the field [8,11] or in the laboratory, on core samples [2,12]. Measurements on cores are more accurate, because it is easier to control the experimental conditions, and a classical method is to perform out-diffusion experiments, whose principle may be described as follows. The core is plunged into a fluid doped with a solute and is left in the solution during a sufficiently long period to evenly saturate the rock. Afterwards, the sample is out-leached by plunging it into a vessel containing the same liquid, but without the solute. Out-diffusion of the solute occurs and the concentration of the external solution is monitored. An interesting variation is to replace the solute by a radioactive element, for instance, a MethylMethAcrylate resin doped in ¹⁴C. With this method, one can obtain autoradiographs of the sample, which basically provide 2D images of the connected porosity. The grey level of each pixel of the image can be converted into a porosity value by comparison with calibrated sources radiographed with the sample (for details, see [3]). An example of autoradiograph and out-diffusion curve from the Kivetty granite (Finland) is shown in Fig. 1. Usually, the out-diffusion curves are interpreted by analytical solutions that are exact in the context of a homogeneous diffusion. In the early times of the experiment, the homogeneous assumption should yield an increasing linear behaviour of the concentration curve with respect to the square root of time [4]. As shown in Fig. 1, this interpretation is not acceptable because the heterogeneous porosity of the rock generates a very complex convex curve.

By using image processing on the autoradiographs, we are able to build calculation grids where each pixel is associated with a known porosity value and an unknown apparent diffusion coefficient. Note that the apparent diffusion is the pore diffusion multiplied by a retardation factor when sorption onto the solid occurs with instantaneous equilibrium. In the following, it is assumed that there is no sorption and, for the sake of simplicity, the pore diffusion will be referred to as diffusion. Numerical simulations of the diffusion can be made, as well as attempts at identifying the values of the diffusion coefficients. By using 2D or reconstructed 3D geometry, one can fit the complex experimental curves and enhance the predictive capabilities of the model, because heterogeneity is accounted for, both in terms of the pore space geometry and the diffusion properties. This work demands cumbersome calculations, because a precise geometry of the pore space requires grids with sometimes several million pixels. The aim of the present work is to propose a new numerical method, much more rapid than classical ones to solve the heterogeneous diffusion. It is also coupled with a mathematical optimisation procedure to find the diffusion coefficients. Finally, a synthetic example is given for illustration purpose.

As mentioned above, the computation grids may count millions of unknowns and to avoid handling the huge and sparse matrices of the Eulerian methods, we have chosen the Lagrangian approach. In this context, the classical Random Walk (RW) method moves particles by jumps in space and over pixels that correspond to the connected porosity [5]. Since the initial concentration of the out-diffusion experiment is constant over the whole porosity, a given number of particles per porous pixel is required to obtain a relevant representation (at least one particle for the less porous pixel of the medium). Moreover, RW needs a minimum of 2–4 jumps for each particle and in each pixel for an accurate representation of the diffusive transport. These constraints argue against the classical RW, leading us to develop a more rapid method using particles in the time domain.

The principle is to calculate the diffusion time needed by a particle to jump from centre to centre of neighbouring pixels. The method is partly inspired by and modified from previous attempts in petroleum engineering to solve the steady-state flow equation [7]. The diffusion equation in heterogeneous media is:

The dimensional analysis of the second term of the right-hand side shows that expression (4) may be rewritten as separate equations, in the form: $\mathrm{d}\mathrm{t}=-({\varphi }_i{V}_i/{\sum }_j{b}_{\mathrm{ij}})\mathrm{d}\left(\log \left({\mathrm{C}}_i\right)\right)$ with an occurrence probability p_ij such as: p_ij=b_ij/∑_jb_ij. Note that b_ij is of dimension [L³T⁻¹], and therefore the probability p_ij corresponds to the fraction of the volumetric flux that passes by diffusion from i to j. In other words, a particle present at the centre of pixel i will jump to the centre of pixel j with a transition probability p_ij and the diffusion time from i to j is given by:

The TDD method is based on a finite-volume formalism and is therefore not free from truncation errors, even if the calculation is performed in time. These errors are of the fourth order with respect to space coordinates or equivalently, of the second order with respect to time. Typically, this error for a jump of size Δx is:

It can be assumed that the image processing of autoradiographs provides the porosity value of each pixel, but the diffusion coefficient remains unknown. Note, however, that the porosity could also be identified, for instance, if the confidence in the image is poor. The inversion might refine the knowledge of the porosity values with the processed image acting as a preconditioning in setting the range of possible values for each grey level. In the following, the problem is simplified by focusing the inversion on the diffusion coefficients. It is obvious that we cannot identify one specific value for each pixel: this would lead to an ill-posed problem with thousands of parameters as compared to limited observations, i.e. the experimental out-diffusion curve. A parameterisation able to strongly reduce the set of the sought parameters is required. In crystalline rocks, the variations of porosity are mostly the consequence of differential weathering and micro-deformations, such as kinks and cracks, which affect the grains, specifically according to their mineralogy. Thus, a strong link exists between the porosity value, the internal pore geometry (constrictivity and tortuosity of pores), and the medium in which it develops. Knowing this, a simple and efficient parameterisation consists in assigning a constant and unknown diffusion coefficient to each grey level of the image. There are as many parameters to identify as grey levels, i.e. a maximum of eight in current applications (the number of thresholds in processing the autoradiographs could be increased, ongoing work shows that 8 to 12 grey levels are sufficient on most images).

The inversion is performed following an indirect approach, i.e. by seeking the parameters iteratively until a criterion of closeness between the observed and the simulated out-diffusion curves is reached. If ${\underset{¯}{d}}^k=\{{\mathrm{d}}_1^k,\cdots ,{\mathrm{d}}_p^k\}$ denotes the vector of the p parameter values at iteration k and ${\underset{¯}{\xi }}^k=\{{\xi }_1^k,\cdots ,{\xi }_n^k\}$ the vector of the n sampled errors between observed and simulated concentrations, the closeness criterion minimises a classical least-square function:

This method has proven its capacity to solve inverse problems with few parameters [6]. The reason may be that the method corresponds to the least-square resolution of the general linear approximation: $\mathrm{f}\left(\mathrm{b}\right)=\mathrm{f}\left(\mathrm{a}\right)+(\mathrm{b}-\mathrm{a})(\partial \mathrm{f}/\partial \mathrm{x})$. Applied to the error vector $\underset{¯}{\xi }$ as a function of the parameter vector $\underset{¯}{d}$, this gives:

When the optimisation is completed, the Gauss–Newton method offers directly, as a by-product, the sensitivity of the model to the parameters. Each term (∂ξ_i/∂d_j) of the Jacobian matrix can be viewed as the local sensitivity (in ξ_i) to the parameter d_j. The terms of J_ξ^T·J_ξ are: ${\sum }_{\lambda =1}^n{\left(\partial {\xi }_{\lambda }\right)}^2/\partial {\mathrm{d}}_i\partial {\mathrm{d}}_j$ and can thus be interpreted as measures of the global sensitivity of the model. The diagonal terms of the matrix are proportional to the variance of the sensitivity to d_i, the out-diagonal terms are proportional to the covariance of the sensitivity to d_i and d_j. It is obvious that the less sensitive the model is to a parameter, the larger the uncertainty on the parameter value. The calculation of the Jacobian matrix is described in the following section.

The example proposed below is a synthetic one, built for illustration purposes and to demonstrate the relevance of the whole inversion procedure, i.e. the TDD transport simulator and the Gauss–Newton method for identifying the parameters. The computation grid is shown in Fig. 2a. It encloses four types of zones with different porosities and diffusion values in a rock matrix, plus the external fluid. Three of these zones are porous patches that could mimic areas of different mineralogy and, as stated before, they are affected by differential weathering (porosity of 0.1, 0.15, 0.25, see Fig. 2). These patches are connected to each other and to the edges of the sample either directly, or by thin stripes that may represent partly filled cracks in the matrix (porosity of 0.3). The total grid size is 100×100 square pixels of 5×10⁻⁴ unit length on a side. In view of the grid, it can be foreseen that the out-diffusion curve is likely to show variations in the slopes of the concentration versus time. This is the logical consequence of the variable connectivity between the pore space and the edges of the sample. Some highly diffusive patches are in direct contact with the external fluid and will contribute most to the concentration curve in early times. On the other hand, the patches with low diffusion or those merely connected to the edges by thin stripes contribute to a delayed concentration increase. The direct problem was solved to get a reference curve for the inversion with the following set of diffusion coefficients: {2×10⁻¹¹,5×10⁻¹³,2×10⁻¹¹,2×10⁻¹²} [L²T⁻¹] ranked according to the indexes of zones in Fig. 2. The initial condition of a constant concentration at any location in the matrix was set up by placing 10 particles in each of the less porous pixels, the other pixels being filled according to this value (e.g., 10 particles for a porosity of 0.1→15 particles for 0.15). With these settings, the direct simulation is performed in less than 10 seconds with a Pentium 2 GHz. The reference out-diffusion curve is shown in Fig. 2b.

To obtain an operational Gauss–Newton inversion, the Jacobian matrix (see Section 3) must be calculated. Contrary to Eulerian models, the Lagrangian framework of the transport simulator (see Section 2) jeopardizes the chance to derive analytically the Jacobian matrix in (12). It was thus processed numerically and the term-by-term evaluation of the Jacobian can be represented by:

In theory, the perturbation should tend toward zero for the best estimation, but the stochastic nature of Lagrangian computations prohibits the use of very small perturbations. With a reasonable number of particles, two direct calculations with the same parameters may show slightly different results. A relative difference of 1% or so on the concentrations is negligible when the problem is to simulate a transport scenario, but for an inversion, it imposes to set the δd_j values so that they induce visible modifications of the concentrations as compared to the uncertainty from the direct calculation. The perturbations were therefore set to 10% of the current parameter values at each iteration. Note that for each perturbed calculation, we did not keep the same chain of generated random numbers u₀₁ (in expression (5)) used to move the particles. The first reason is that a perturbation on a parameter modifies both the transition probability and the transition time distribution from pixel to pixel. Thus, keeping the same random number for a transition has no physical consistency. In addition, this procedure has shown to be cumbersome in terms of storage capacity and there is no significant gain in computation costs.

The starting point of the inversion was set at 5×10⁻¹² [L²T⁻¹] for any porous zone and the diffusion coefficient for the external fluid was assumed to be known (10⁻¹⁰ L²T⁻¹). The Gauss–Newton algorithm finds its optimal solution in 12 iterations and yields the following diffusion coefficients {1.85×10⁻¹¹,1.62×10⁻¹²,1.52×10⁻¹¹,3.45×10⁻¹²}. Whereas areas 1, 3, and 4 are found with precision, area 2 is poorly identified. The reason is that the particles from area 2 contribute weakly to the out-diffusion curve, since (1) the area is not well connected to the edges of the sample and initially encloses few particles, (2) the low diffusion coefficient causes these particles to be present in the flat part of the concentration curve (latest times) where the model is not sensitive. This is confirmed by the diagonal of the Jacobian matrix, whose values have been normalized by the optimal parameters to produce something with the dimension of a quadratic error: (∑_λ=1ⁿ(δξ_λ/δd_i)²)d_i². These values are {6.16,2.88,13.4,45.2}×10⁴ [ML⁻³]² and show the model to be less sensitive to parameter #2. Making the starting point vary does not really change anything, except when it is very close to the true solution. This feature is inherent in the inverse problem and not the consequence of the Gauss–Newton optimiser. However, some out-diffusion curves may show different portions, each of them stemming from particles of fairly well identified areas. By adding a flag to these particles, one could focus the inversion on the parameter that is assumed to mainly influence a specific part of the curve. This would increase the speed of convergence and probably, the sensitivity to parameters. We also carried out inversions in a more classical way using an Eulerian finite volume formalism, a Gauss–Newton optimiser and an analytical derivation of the Jacobian matrix (which is straightforward in the Eulerian context). The performance in terms of parameter identification is similar to TDD inversion. However, each direct simulation requires ten times as much computation time as TDD. Even though the analytic Jacobian lowers convergence iterations from 12 to 8, the TDD inversion remains much more rapid than a classical discrete method.

Out-diffusion experiments are widely used to assess the safety of potential underground repository sites. We have shown that a comprehensive interpretation of the diffusion through a heterogeneous rock matrix needs intensive computations. We built a Lagrangian simulator in the time domain that is rapid and well suited to transport simulations over grids extracted from autoradiographs of rock samples. The calculations are well conditioned to both the spatial distribution and the local values of the porosity. In addition, the coupling with a Gauss–Newton optimiser makes it possible to seek automatically the diffusion coefficients and lowers the subjectivity of a classical fitting by ‘trial and error’. Further works will focus on how to increase the sensitivity of the inversion to the parameters. For instance, some theoretical developments are in progress to derive the analytical form of the Jacobian matrix and calculate its coefficients by a method similar to TDD. Basically, this improvement does not modify the principle of the Gauss–Newton inversion, but a better estimation of the descent direction would probably lower the number of convergence iterations. Other lines of research also deserve some interest, for instance, adding mechanisms such as sorption with first-order kinetics and handled in the time domain.