Probabilistic modelling of calcium leaching in a tunnel for nuclear waste disposal

This works aims at studying, through the Monte-Carlo method, the influence of the leaching kinetics parameters spatial variability on the lilfespan of a concrete structure. The considered structure is a tunnel for nuclear waste disposal. It is observed that the expected value for the lifespan estimated when considering the material spatial variability is significantly lower than the lifespan estimated with a single simulation considering uniform parameters fields, equal to the expected value for each parameter.


INTRODUCTION
Deep nuclear waste disposal facilities need to be studied over periods of one or more order of magnitude greater than those of classical civil engineering.This leads to consider some degradation phenomena that are most of the time not taken into account because of their slow kinetics.This is the case with calcium leaching for concrete structures exposed to water from the host rock in a deep underground disposal.Moreover, civil engineering constructions are naturally subjected to variability, from many origins: their dimensions, their construction processes, their exposure to several loadings...Among these variables, which influence the behaviour of the structure and, thus, its lifespan, some of them introduce a variability on the material characteristics, evolving with time and not necessary homogeneous in the whole structure, whereas some others create a variability on the loading itself, which the structure is subjected to.
The structure which we study is a tunnel for radioactive waste disposal.The objective is to confine the waste in a deep geological formation to protect the environment from the radionuclides they contain.This confinement has to be guaranteed over large time scales (up to several hundreds of thousands years) without requiring maintenance.The waste is contained in packages of metal or concrete, depending on the nature of the waste.These packages are arranged in cells, being disposal cavities in the rock.Under the hypothesis of reversible storage, packages are simply stored in the cells.If the option of underground disposal is finally adopted, the cells, and the tunnels leading to them, will be sealed.
The structure we consider is one of these cells.For simplicity reasons and to ease the calculations, the considered geometry for the structure is very schematic: we simply consider a hollow cylinder, with an external radius of 8 m, and a thickness of 1 m.The disposal cell can reach a length of several hundreds meters.Figure 1 schematically represents the problem that we model.A tunnel for the waste packages disposal is dug in the soil, and a concrete wall protects the environment from the waste containers contained within the tunnel.The site water contained in the geological formation exposes the concrete tunnel to the risk of leaching.Even though this kind of degradation is very slow, it is to be considered as the required lifespan for such structures is several hundreds of thousands years.A numerical model for calcium leaching is developed and implemented in Finite Volumes simulation code.This model will be briefly presented in the first section of this paper.The second section introduces the statistical data we use to model the variability of our input parameters to model calcium leaching, being water porosity and coefficient of tortuosity.The last part is devoted to the investigation about the influence of this spatial variability on the lifespan of the tunnel.This study is performed by carrying out the Monte-Carlo method, which requires a large number of calculations.Moreover, as the accuracy of estimating the position of the degradation front depends on the refinement of the Finite Volumes mesh, it has been decided not to model the entire cell disposal, but only a part of it.The section considered in this application corresponds to an angular fraction of /12, as it appears on Fig. 1.This is close to a block of 1 m high, wide and long.
Let us specify here that this problem is treated under a few strong assumptions, besides the geometric simplification of the structure: • the influence of temperature is neglected, that is to say that the simulations are performed at constant temperature throughout the duration of the degradation (this assumption leads to overestimate the lifespan of the structure since the temperature might reach 60 • C for 10 years and 50 • C for 100 years according to the type of waste); • any other interaction with the environment than leaching in pure water is not taken into account in modelling; • leaching degradation is supposed to occur in pure water, whereas site water is actually loaded with different ions; • no coupling with the mechanical behaviour of the structure is taken into account.Finally, let us specify that the so-called lifespan of the structure in this study is the time necessary for the portlandite dissolution front to reach the inner wall of the tunnel.

CALCIUM LEACHING MODELLING
The phenomenon of leaching consists in the dissolution of solid calcium in cement hydrates when concrete is exposed to any aggressive solution (most of the time pure water or at least water with low calcium concentration).Indeed, concrete porous solution is very basic (pH around 13) and several ionic species are highly concentrated [1].Therefore, water with low mineral concentration is an aggressive environment for concrete.Concentration gradients between external environment and porous solution induce diffusion of the main ionic species (alkalis, calcium, hydroxides) from the porous solution to the external environment.Consequently, the initial chemical equilibrium between the porous solution and 04004-p.2AMP 2010 the solid phases of the material is modified.A new local equilibrium is then reached with the dissolution of solid phases [2].
Considering that the dissolution is instantaneous (local equilibrium) and that only calcium species are to be taken into account, the leaching of a cement paste can be described by the mass balance equation of calcium (1), as it was proposed by [3], where S Ca is the solid calcium concentration, C Ca is the liquid calcium concentration, D is the calcium effective diffusivity in porous material and is the porosity.One can recognise the two main phenomena involved in the leaching process: on the one hand, the solid and liquid phases of calcium in the cement paste and in the porous solution are in chemical equilibrium; on the other hand, the ionic calcium species diffuse through the material porosity.
The non-linearity of the former equation is mainly due to the diffusivity which depends on the porosity, itself depending on the solid calcium concentration, and on the non-linearity between S Ca and C Ca .The Finite Volume method [4] seems to be particularly accurate for our concern.Indeed it is well adapted to diffusion problems because its equations remain conservative.This is why a Finite Volume scheme has been chosen to pursue this study.
Considering the chemical laws of local equilibrium between the concentrations of liquid and solid calcium, every quantity appearing in (1) can be derived from the concentration of solid calcium, which reduces the problem to a single variable S Ca .Considering concrete rather than cement paste leads to the introduction of a tortuosity parameter standing for the possible opposite effects of aggregates: on the one hand aggregates are areas where no diffusion occurs, but on the other hand they introduce an Interfacial Transition Zone where the diffusivity can increase due to a larger porosity of the ITZ in comparison with bulk cement paste [5].This approach was successfully used by [6] and appears in (2) in order to express the diffusivity of calcium species through concrete depending on the cement paste porosity.
The implementation and validation of this numerical model are presented in [7]; this latter study also enhances that the most influential parameters on the leaching phenomenon kinetics are the porosity and the coefficient of tortuosity.
The main issue here is that porosity is a measurable parameter which can directly be observed for any considered concrete sample whereas the coefficient of tortuosity is a non-physical parameter and depends on the numerical model we chose.Therefore, an inverse analysis tool has to be developed so as to identify this coefficient from the observed experimental data, being the degradation depths at several terms of an accelerated leaching test.Among the numerous tools developed for inverse analysis, it has been chosen to use an Artificial Neural Network (ANN).It is used as follows: the input data are the values measured for the porosity and the degradation depths for each sample, and a mean value of temperature for each period of degradation.The output data is an identified value of coefficient of tortuosity.

MONTE-CARLO INTEGRATION AND RANDOM FIELDS GENERATION
Dealing with statistics computations for some outputs of interest, the Monte-Carlo fundamental principle is to estimate the integral -i.e. the expected value -of this functional through its evaluation for a finite number of realisations which are randomly chosen [8].It is worth noting that despite its rather slow convergence, the Monte-Carlo method may be as accurate as required provided the number of "integration points" is large enough.For every realisation, the method consists in solving a deterministic problem through the Finite Volumes model presented in Section 2.

04004-p.3 EPJ Web of Conferences
As explained above, the key step to carry out the Monte-Carlo integration is to generate an important number of realisations of the porosity and the coefficient of tortuosity random fields.Although these realisations must be independent, they are all spatially correlated and the synthesis process is not trivial.A very efficient way to achieve it is to employ the Karhunen-Loève expansion as explained by [9].The key idea of the Karhunen-Loève expansion is to expand any second-order random field f ( x, ) into a serie (3).Each term of this serie is a product of two terms, the first one depending on the spatial variable x and the second one depending on the stochastic variable .Thus, the main feature of such a decomposition is the separation of variable. ( The first term of ( 3) is the expected value of the random field.The stochastic dependency appears through the random variables i, while the spatial dependency is introduced through the covariance kernel eigenmodes ( i, i ), where the i are the eigenvalues and the i are the corresponding eigenvectors.Moreover, it is worth noting that the random variables i form a set of orthonormal variables.Finally, dealing with Gaussian random fields, the i are Gaussian random variables as well and, being uncorrelated, they are independent.Thus (3) provides a very efficient way to synthesise correlated random fields realisations, independent from one another.
Considering the generalised eigenvalues problem, we are interested in determining the largest eigenvalues.A convenient way to achieve this is to perform a Finite Elements discretisation of the problem (see [10]), leading to the discrete generalised eigenvalues problem ( 4), where M is the mass matrix with a unitary mass density and C is the covariance matrix.
This numerical problem is very alike to eigenvalues problems encountered in structural dynamics.Therefore, a Finite Element code (FEAP) has been used to determine the Karhunen-Loève eigenmodes by re-using the FEAP routines with the appropriate matrixes.It is worth noting that, even though M is a sparse matrix -due to the FE properties -the covariance matrix C is a full one.This point may lead to some difficulties linked to disposal and memory requirements.Once a set of the largest eigenvalues, as well as the corresponding eigenvectors, has been determined, the truncated Karhunen-Loève expansion enables the approximation (5) for any second-order random field.In this approximation m is the number of modes which are considered and it can be shown that ( 5) is optimal along m [11,12].Assuming the random field to be Gaussian, the random variables i are normally distributed with zero expectation and variance equal to one.This series truncation introduces an error in the field approximation, which decreases when m increases.The theoretical background of the Karhunen-Loève expansion can be carried out to calculate the eigenmodes of the covariance kernel.Once these modes are in hands, the expansion can be exploited as many times as necessary to compute the random fields realisations required for the Monte-Carlo Method.

VARIABILITY OF POROSITY AND COEFFICIENT OF TORTUOSITY
The work presented here is part of a project named APPLET and funded by the National Research Agency (France).This project aims at investigating the variability of the characteristics of concrete, such as its mechanical properties or durability indicators [13].The main purpose is to acquire statistical data on concrete properties through a large experimental campaign.This campaign is led in partnership with the Vinci Company, and consists in following two real construction operations, corresponding to 04004-p.4 AMP 2010 two different kinds of concrete mix design, the first one being a high-performance concrete with fly ash and the second one being an ordinary concrete.
For each construction site and so for each concrete formulation, 40 batches are characterised through several tests performed in laboratory (such as compressive and splitting tensile strength, static and dynamic Young modulus, electrical resistivity, etc.).All the concrete samples necessary for the tests are directly provided by the construction site to the involved laboratories.The concrete samples are protected from drying thanks to plastic bags.Then they are stored in saturated lime water.The concrete specimens are one year old at the beginning of each test.This campaign aims at acquiring statistical data on the material characteristics and investigating the variability of these characteristics within several batches of a same formulation and also between different concrete mix design.The porosity measurement and the accelerated leaching tests were performed in LMT Cachan.The whole experimental campaign is presented in detail in [14].For this study, only the results on the A1 concrete are considered, because the following numerical investigation is based on the experimental observations for this concrete mix design.
The so-called A1 concrete is the high-performance concrete with fly ash used for the A84 highway tunnel near Paris (France).Table 1 presents the concrete mix design used for this building operation.In Table 2, the basic physical characteristics of this material are summarised.In this table appear the average values experimentally observed during the APPLET campaign for the electrical resistivity el , the volumic weigth , the compressive strength 28 days after casting F c 28 and the compressive strength after one year in saturated lime solution storage F c 365 , the tensile splitting strength Fc and the Young modulus E.  As mentioned above, one of the main parameters involved in the leaching process is the water porosity measured on sound concrete.Each batch of concrete studied in the project is characterised by the measurement of its water porosity, which is calculated on the basis of three mass measurements: two performed on a saturated slice of concrete (full atmospheric and hydrostatic weight, M air and M water ) and one performed on the dried slice (105 • C in a furnace till constant mass, M dry ) [15], according to equation (6). 04004-p.5

EPJ Web of Conferences
For each batch, a concrete sample is devoted to a leaching test.Leaching in pure water usually presents a very slow kinetic and one of the most common ways to accelerate it is to use ammonium nitrate solution [17,18].Calcium concentration increases significantly (from 22 mol/m 3 to 2700 mol/m 3 for a concentration of 6 mol/L of ammonium nitrate), which accelerates the degradation kinetics at least by a factor 100.At fixed dates (28, 56, 98 and 210 days), the degradation depth is revealed with phenolphthalein.
Actually, it seems that the degradation depth revealed with phenolphthalein is not exactly the position of the portlandite dissolution front [16], but the ratio between both is not completely acknowledged: it does not seem constant in time neither with different degradation conditions.Therefore, it will be considered in this study that the degradation depth is the one revealed with phenolphthalein.
For each concrete sample, the coefficient of tortuosity is identified with the Artificial Neural Network introduced in Section 2. Thanks to this experimental and identification campaign, one can quantify the variability observed between several batches of the same concrete mix design on a building site for porosity and coefficient of tortuosity.These results are presented for the A1 concrete of the APPLET project in Table 3.
Table 3. Water porosity and coefficient of tortuosity: number of tested samples, mean value and coefficient of variation for A1 concrete of the APPLET project.

Nbr. Mean value CoV [%] Mean value CoV [%] A1
40 12.9 7.9 0.134 15.1 If the experimental campaign results presented above allow us to propose a quantification of the variability for both porosity and coefficient of tortuosity (especially in terms of average and standard deviation), they do not, however, provide any information about the spatial correlation for these fields.Therefore, an experimental wall was cast, and samples were taken along two horizontal lines and one horizontal line.This wall was made with a construction joint at half height, so that each horizontal line is in a different batch.The samples are cores of 10 cm large and 30 cm high.However, the two batches were cast the same day, so there is very few difference between the two batches and the experimental results shall be regarded as from one batch only.On each sample, water porosity and resistance to an accelerated leaching test are measured, following the same protocol as presented above.So we can identify, for each sample, the water porosity and the coefficient of tortuosity .
In our case, we assume that our random fields (porosity and coefficient of tortuosity) are spatially correlated according to the covariance function given in equation (7) where C(x, y) is the covariance between two points, x and y, and 2 is the random field variance, while L c is a scalar parameter called correlation length.This latter parameter reflects the importance of the field spatial correlation: the higher the value, the stronger the field is correlated, so that for an infinite correlation length the field would tend to be uniform.
This covariance function has the advantage of being differentiable at 0. For better understanding, let us recall that the covariance is an indicator of the degree of dependence between two variables: if the covariance of two variables is zero, they are independent from one another, that is to say that their respective variations are independent.Thus the parameter L c translates the random fields spatial correlation, and the experiments on the experimental wall aim at identifying this correlation length.Provided one has several realisations of a random field, it is possible to calculate the covariance matrix on the considered geometry, meaning to calculate the covariance between each pair of points of the 04004-p.6 geometry and to identify the correlation length.Indeed, by definition, the covariance between two points x and y for a random variable respectively equal to X and Y , written as where the expected value of a random variable is denoted E. However, if one has only a single realisation of the considered random field (as it is the case here: one experimental wall), it is possible to determine the covariance from a variogram.For a given random variable, the value of the variogram at distance d is the half-mean square differences of its realisations on the set of points spaced of d.For a given random field, if covariance is defined, its value for a distance d is equal to the difference between the covariance at 0 and the value of the variogram at d. Covariance at 0 simply is the variance.Since the variogram can be determined from a single realisation of the field, discrete values of the covariance can be deduced.Then we identify the value of the parameter L c in the covariance function (7).
To conclude, Fig. 2 shows the covariance values obtained from the measurements on the wall compared to the corresponding covariance function.One can plot the normalised covariance (with regards to the variance for each variable: porosity and coefficient of tortuosity) on the horizontal line in batch 2 for instance (see Fig. 1 Fig. 2).
It is observed in all cases that the correlation lengths are by the order of magnitude of one meter, and are similar for all three studied lines, whatever the considered variable.However, since the accordance between the proposed correlation function and the variogram values observed for the two variables is only relatively satisfactory, and given the low number of measurement points from which the study was conducted and the small size of the experimental wall (with regards to the scale of an entire structure), these values of correlation lengths (summarized in Table.4) should be considered with caution.

RESULTS AND ANALYSIS
The imposed input spatial variability is observed on the samples from the experimental campaign (cf.Section 4).   .These two variables are assumed independent.Three values of correlation lengths are tested, by the order of magnitude identified on the experimental wall: 0.5 m, 1 m and 2 m.
The lifespans estimated by taking into account the input parameters variability are to be considered with respect to a reference value which is given by a deterministic simulation, that is to say, a simulation with uniform fields of porosity and coefficient of tortuosity, equal to the expected value of these parameters.Indeed, the expected value for any output of a problem is not equal to the output of the same problem considered with the expected values of its input parameters.Figure 3 shows the mesh chosen for these simulations.The sides of the element, which are not exposed to site water, are assumed to be subjected to a null flux.There are 10 control volumes along the length of the structure (parallel to a generatrix of the cylinder), 10 orthoradial volumes and 20 volumes in the thickness.
The reference simulation is conducted with uniform fields of porosity and coefficient of tortuosity, being respectively equal to 12.9% and 0.134.The result in this case is the reference value for the lifespan of the tunnel: 254,000 years Let us observe the results in the case of a correlation length of 0.5 m. Figure 4 shows the results of the Monte Carlo simulations and compares them with the result of the reference simulation (Figure 4(a)).Figure 1 represents the lifespans distribution and it is to be noticed that the expected value for the probabilistic lifespan is significantly lower than the reference value (227,554 years vs. 254,000 years).In this figure, one can also observe a generalised extreme value distribution law, which is the statistical law which reproduces better the results distribution.structures have a lifetime of less than 198,760 years, to be compared to the 254,000 years of the reference simulation, which stands for a difference of almost 22%.
It is verified that in all cases of correlation length, the expected value for simulations taking into account the input fields variability is significantly lower than the reference value of lifespan: 245,648 years and 225,883 years respectively for a correlation length of 1 m and 2 m, instead of 254,000 years for uniform fields.
It also is noticeable that the coefficient of variation increases with the correlation length: 10.6% for L c = 50 cm, 12.0% for L c = 1 m and 16.4% for L c = 2 m.This increase in the coefficient of variation with the correlation length is mainly due to the fact that a large correlation length (in the latter case the correlation length is two times greater than the tunnel thickness) implies the existence of "extreme cases" among the realisations of the random field, that is to say uniformly low (or on the contrary high) realisations, which would lead to penetration velocity of the degradation front substantially smaller (or larger) compared to the real uniform case with the expected values of porosity and coefficient of tortuosity.
For each case of correlation lengths, 3 distribution laws are fitted on the results distribution: lognormal, logistic and generalised extreme value distribution law.All three laws reveal to be quite relevant to model the results distribution.The values of the statistical laws parameters, identified on the simulation results for the three studied cases of correlation lengths, are listed in Table 5. Figure 5 compares the simulation results distribution for different values of correlation length with the three distribution models (log-normal, logistic and generalised extreme value).Note that from a safety point of view, it is important to choose the most relevant distribution model for the shortest lifespans.Thus the most accurate models are the generalised extreme value distribution for the correlation lengths of 50 cm and 2 m, and the logistic distribution for the correlation length of 1 m.Also note that in all cases, the reference value of the lifespan is much higher than the median 04004-p.9

Figure 1 .
Figure 1.Schematic representation of the considered structure: part of a tunnel for nuclear waste disposal exposed to ground water.

AMP 2010 Figure 2 .
Figure 2. Correlation lengths identification from the variograms for the two considered random variables on the horizontal line in batch 2.

Figure 3 .
Figure 3. Finite Volumes mesh for the modelled part of the tunnel.

Figure 4 .
Figure 4. Fitting of a generalised extreme value distribution law parameters in the case L c = 0.5 m.

Figure 4 (Figure 5 .
Figure 5.Comparison between lifespans cumulative distribution and log-normal, logistic and generalised extreme value cumulative distribution laws.

Table 1 .
Concrete mix design for A1 building operation.

Table 2 .
Basic physical characteristics observed for A1 concrete.

Table 4 .
Identified correlation lengths for porosity and coefficient of tortuosity on the experimental wall.
Table 3 identifies the variability observed for porosity and the coefficient of tortuosity

Table 5 .
Identified values of the distribution laws parameters.