Simulation of Reactive Diffusion in Clays By a Continuous-Time

Furthermore, the Markovian nature of the modeling approach lends itself to a continuous-time particle-tracking scheme of resolution of the model equat...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCC

Simulation of Reactive Diffusion in Clays By a Continuous-Time Markovian Particle-Tracking Scheme Francesco Cadini*,† and Enrico Zio†,‡ †

Dipartimento di Energia, Politecnico di Milano, Via Ponzio 34/3, I-20133 Milan, Italy Chair on Systems Science and the Energetic Challenge, European Foundation for New Energy, Électricité de France Ecole Centrale Paris and Supelec, Grande Voie des Vignes 92295, Chatenay-Malabry Cedex, France

‡

ABSTRACT: Clay minerals and clay rocks are considered as efficient components of engineered and natural barriers in many high-level radioactive waste disposal programs worldwide because of their low permeabilities and high sorption capabilities. In this paper we present an approach for modeling solute diffusive transport in saturated clay minerals at the structure map scale based on an extension of the Kolmogorov−Dmitriev theory of stochastic branching processes. The proposed modeling framework allows a simple description of (i) the dual-porosity behavior typical of some compacted clay structures, where diffusion occurs both in the interlayers and in the pore space between clay particles, and (ii) the chemical interactions between the dissolved ions and the solid matrix surfaces (both at nonequilibrium and equilibrium conditions). Furthermore, the Markovian nature of the modeling approach lends itself to a continuous-time particle-tracking scheme of resolution of the model equations, which, in general, easily allows us to deal with the complex geometrical structures of the porous media. Finally, in order to account for the uncertainty in the geometrical properties of the clay samples within the proposed modeling framework, the structure map is represented as a stochastic network of interconnected, one-dimensional interlayers with different lengths and orientations. The combined methods are demonstrated on a two-dimensional literature case study of diffusion through a compacted clay sample, and the resulting average diffusion coefficients are in general agreement with those obtained by a different simulation-based technique found in the literature. Then, the methods are applied to model the diffusion of the weakly sorbing cation Na+, providing results in accordance with physical intuition.

1. INTRODUCTION Clays and clay rocks have been and still are under investigation as major components to be used in toxic waste deposits and deep geological radioactive waste repositories.1−4 In particular, compacted clays, such as bentonite, are widely considered to be the optimal materials for the engineered (or artificial) barriers, whereas natural clay rocks, i.e. fine-grained, porous materials containing a significant proportion of clay minerals,2 are envisaged in many countries as potential host rock formations for disposal facilities.3,5−7 The interest in these materials is mainly due to their extremely low permeabilities and large sorption capabilities of cationic species by the negatively charged clay structures.2−4,8 Bentonite is essentially composed of montmorillonite, which is a smectite-type clay mineral, whereas clay minerals in clay rocks are in general made up of both illite- and smectite-type minerals in different amounts. Both minerals are made up of platelets of positive charge deficient alumino-silicate layers stacked one above the other. The resulting net negative charge is balanced by cations located in the interlayer pore space and attached to the external aggregate surfaces. In smectites, differently from illites, interlayer cations and their associated hydration water can exchange with cations and water in the pore solution,2 providing these minerals with the capability of © 2013 American Chemical Society

swelling in the presence of water. As a consequence, smectitetype minerals exhibit in general extremely low permeabilities, and molecular diffusion can be assumed as the dominant mass transport process.9 The diffusion of water or solute particles in smectite-type minerals is thus strongly influenced by the geometric layout of the complex microstructure of these clays, described in Glaus et al.9 as a dense “house of cards” arrangement of clay platelets, which are in turn composed of stacks of clay layers. At the same time, water can be found in two different physicochemical states in smectites: the free water in the pore space between the edges of the clay platelets, whose properties are similar to those of bulk water, and the interlayer water present between the clay layers, whose properties are influenced by the counterbalancing cations.9 Thus, in general these materials behave as doubleporosity media, in which water particles can migrate in the macropores filled with free water and the nanopores of the interlayer water, each of which has different diffusive properties; as the level of compaction increases, however, macropores tend to reduce their volume and interlayer diffusion seems to Received: June 4, 2013 Revised: August 5, 2013 Published: August 7, 2013 18510

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

become the dominant pathway.9 It is worth noting that the high density of negative surface charge, especially in compacted clays, allows diffusion of only neutral and cationic species in the interlayer water, whereas anions are excluded because of the strong electrostatic repulsion of the clay layer surfaces.2,9 Finally, the diffusion of cations is further influenced by chemical interactions with the surfaces of the solid phase, generally due to the combined effect of ion exchange and surface complexation reactions.8,9 Under growing international pressure on the definitive disposal of high-level radioactive wastes, much effort has been placed on improving the understanding of diffusion-driven transport in clays and clayrocks, from both the experimental and the analytical/numerical points of view: the interested reader may refer to Altmann et al.2 and Robinet et al.8 for thorough reviews of recent works and the state of the art with a focus on the issue of cation diffusion. Despite the large availability of diffusion data and models, many aspects are still under investigation, especially with regard to (i) the still poor understanding of the phenomena involved in the diffusion process of charged particles (cation surface diffusion, for example),3,9 (ii) the consequent need to improve the numerical models in order to account for these phenomena within the complex geometrical structures of clays, and (iii) the up-scaling of local scale numerical models for the prediction of mass transport on larger scales and longer times. Recently, Churakov and Gimmi4 presented a numerical modeling approach for the diffusion of water or ions in clays at the structure map, or even larger, scale. The method stems from the up-scaling of molecular dynamics simulations of diffusion at the molecular scale to classical, continuous-space, and discrete-time random walk simulations10−12 at the larger scale. The up-scaling procedure allowed the capture of the details of the surface molecule interactions by the molecular dynamics simulations and, at the same time, accounted for the complex geometries of the clay minerals by the random walk simulations. However, only the mobility in porous media of water or of aqueous ions at equilibrium was considered, without taking into account the dynamic evolution of a nonequilibrium system toward equilibrium due to strong site-specific sorption interactions.4 In order to tackle this issue, in this paper we present a particle-tracking-based approach, an alternative to the classical random walk method, for modeling solute transport in clay minerals at the pore structure map scale. The approach is based on an original extension of the Kolmogorov−Dmitriev (KD) stochastic branching model13 applied to the transport of contaminants in groundwater.14 The model aims at describing the stochastic transport of the populations of solute particles in the clay pore spaces by means of a discrete-state, continuoustime Markov process based on transition probability distribution functions defined on the basis of physical reasonings. Partial differential equations, known as the forward Kolmogorov equations (FKEs), are then derived for the expected values of the different particle populations, which are proportional to the solute concentrations. The FKEs are then solved by a particle-tracking (PT) algorithm which can be simply derived from the underlying Markovian formalism. The resulting Kolmogorov−Dmitriev-based particle-tracking (KDPT) algorithm allows for straightforward accounting for the complex geometric structure of the clay particles, as in the classical random walk method, and at the same time also for the

nonequilibrium dynamics of the chemical interactions between the dissolved species and the surface. The capabilities of the KDPT approach are shown with reference to a case study of diffusion within compacted clay samples made up of montmorillonite or pyrophillite taken from Churakov and Gimmi.4 Differently from Churakov and Gimmi,4 we present a method for sampling, within each particle-tracking simulation, realizations of the two-dimensional structure map by building a “network” of connected onedimensional domains with stochastic lengths and orientations, representing both the interlayers and the pore spaces between edges. Thus, geometrical uncertainties are automatically taken into account in the estimation of the average diffusion coefficients. The paper is organized as follows. In Section 2, the KDPT method is briefly recalled and customized to the case of diffusive transport of a solute in the pore space of clay minerals in the presence of chemical reactions with the solid surface. In Section 3, the performance of the KDPT method is demonstrated on the random network of interconnected onedimensional interlayers built on the basis of the information available in Churakov and Gimmi.4 Conclusions on the capabilities and limitations of the proposed modeling scheme are drawn in Section 4.

2. KOLMOGOROV−DMITRIEV PARTICLE-TRACKING MODEL In this section, we briefly recall the continuous-time PT algorithm based on the KD model of contaminant transport in a porous medium (KDPT).14−17 The presentation is tailored to the case of diffusive transport of solutes in clays. The key feature of KD modeling is that different types of particles, characterized by different stochastic behaviors, are introduced to represent the solute in its possible different states (physical or chemical) and regions of space (positions). Partial differential equations, known as the FKEs, are then derived for the expected values of the different particle populations, which are proportional to the solute concentrations. Within a Markovian description of the stochastic space−time evolution of a system of different particles, the FKEs of the KD model are naturally suited to a PT-based solution. The resulting KDPT algorithm amounts to simulating a large number of realizations of the migration fates of contaminant particles by sampling from the probability density functions of the model. The KDPT benefits from the well-known advantages of Lagrangian methods over Eulerian methods, i.e, (i) avoiding the computational burden arising from the need to handle huge and sparse matrices, especially when dealing with strongly heterogeneous porous media18 and (ii) automatically preserving the total solute mass. In addition, by the introduction of new particle types, the discrete-state nature of the KDPT method allows a straightforward and simultaneous inclusion of many physical and chemical processes, as shown for example in Marseguerra and Zio14 and Cadini et al.15 For this, no additional modeling efforts are required, but additional model parameters are introduced. Finally, the KDPT is continuous in time; thus, adaptive time stepping, which is required by spacebased methods, is avoided and a larger computational efficiency in the presence of heterogeneous media is potentially achieved.19 In what follows, the main features of the KD model will be recalled with reference to a one-dimensional domain,20 although the method can be easily extended to deal with 18511

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

solute transport to identify the transition rates of the KD model. The one-dimensional diffusion-reaction equation for the solute concentration reads4,21

two- or three-dimensional problems. The domain is subdivided into Nz discrete zones, z = 1, 2, ..., Nz. Several categories of particles can be introduced according to the different states in which the solute particle under analysis can be found. Thus, the system is in general made up of m = Nc· Nz different kinds of particles, representing the Nc solute particle types at their location in one of the zones z = 1, 2, ..., Nz. Each particle may disappear or give rise to one particle of the m − 1 remaining kinds according to given transition probability laws. It is assumed that (i) the stochastic process is Markovian, i.e., a particle of the kth kind, k = 1, ..., m, gives rise to a branching process independently of its past history; (ii) the process is linear, i.e., the particles do not interact with each other; and (iii) within a generic time interval (t, t + dt), with dt sufficiently small, only one particle transition may occur. Within this framework, in our case we consider Nc = 2 particle types: the solutons, representing the solute particles migrating within the saturated clay pores, and the trappons, representing the solute particles sorbed on the clay surfaces. The mass transport process is then mimicked as follows. In the time interval dt, a soluton can undergo a sorption transition in zone z, thus transforming into a trappon in zone z; the probability of occurrence of this transition, conditioned on the fact that the particle was a soluton at time t in zone z, is ls→t,z(t)dt, where ls→t,z(t) is the corresponding transition rate. Alternatively, within the same time interval dt, the soluton can travel to one of the neighboring zones, z + 1 or z − 1, with transition rates fs,z(t) (forward) and bs,z(t) (backward), respectively. On the other hand, a trappon in zone z can be released into the solution and become a soluton in zone z with transition rate lt→s,z(t). For simplicity and with no loss of generality, we now further assume that (i) the underlying Markov process is homogeneous in time, i.e., the probability density functions governing the stochastic behavior of the system do not change in time and (ii) the medium is homogeneous; thus, all the transient rates become constant in time and space, i.e., ls→t,z(t) → ls→t, lt→s,z(t) → lt→s, fs,z(t) → fs, and bs,z(t) → bs. To describe the space and time evolution of the system of particles, it is possible to write the following system of 2Nz coupled differential equations (forward Kolmogorov equations) for the evolution of the expected values of the numbers of solutons (S) and trappons (T) at time t in zone z′:14 ∂S(z′, t |1s , z , 0) ∂t

∂C ∂ 2C = Dil 2 − f (C , F ) ∂t ∂x

where C is the solute concentration in water (mass of solute per a generic volume of water), F the concentration of solute sorbed on the solid matrix interfaces (mass of solute per mass of solid), Dil the microscopic diffusion coefficient in the interlayers, and f(C,F) the volumetric rate of exchange of solute mass between the liquid and the solid phases (isotherm), which is in general a nonlinear function of C and F. Equation 2.a must be coupled with a mass balance equation for F f (C , F ) ∂F =− ∂t ρs

(C(z + 1, t ) + C(z − 1, t ) − 2C(z , t )) ∂C = Dil ∂t Δz 2 + f (C(z , t ), F(z , t )) f (C(z , t ), F(z , t )) ∂F =− ∂t ρs

∂t

fs = bs =

(3.b)

Dil Δz 2

(4)

f (C(z , t ), F(z , t )) = ρs [lt → sF(z′, t |1s , z , 0) − ls → tS(z′, t |1s , z , 0)]

(5)

The identification of the values of the exchange rates ls→t and lt→s lies beyond the scope of this work, but their effect on the model outcomes will be discussed in the next section. Given proper boundary and initial conditions, the solution of eq 1.a, with the parameters of eq 4 and for fixed values of the exchange rates ls→t and lt→s, can be found by resorting to the continuous-time particle-tracking scheme KDPT described in Cadini et al.15 The method is derived from the underlying Markovian description of the stochastic evolution of the system of particles, whereby the stochastic migration of a large number M of solute particles is simulated by repeatedly sampling their births from the release sources and their transitions across the medium compartments. During each simulation, the positions of the particle are tracked so that, at the end of the M simulations, the corresponding expected values can be estimated by computing the sample means of the positions recorded at the predefined time instants. This procedure will be further detailed in the next section with reference to the simulation of solute diffusion in clays.

(1.a)

= −lt → sT (z′, t |1s , z , 0)

+ ls → tS(z′, t |1s , z , 0)

(3.a)

These equations are formally identical to eqs 1.a and 1.b, and by a term-by-term comparison it is possible to write the following relationships for the forward and backward rates of the solutons and for the volumetric rate of exchange of solute mass between the liquid and the solid phases:

+ bsS(z′ + 1, t |1s , z , 0) + fs S(z′ − 1|1s , z , 0)

∂T (z′, t |1s , z , 0)

(2.b)

where ρs is the solid phase density or bulk dry density. Upon discretization of the above equations in space by a centered Euler finite difference method, we get

= −[bs + fs + ls → t , z ′]S(z′, t |1s , z 0)

+ lt → sT (z′, t |1s , z , 0)

(2.a)

(1.b)

with z′ = 1, 2, ..., Nz, and where the initial condition of one soluton in zone z at time t = 0 (1s,z) is assumed. For any timeand space-dependent contaminant source, the solution of the above system of equations can be found by convolution, thanks to the process linearity. In analogy to the approach proposed in Ferrara et al.20 and Cadini et al.,15 a term-by-term relation can be set up with the classical mass transport equations of one-dimensional diffusive 18512

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

Figure 1. Mean squared displacements of a tracer and associated ±1σ uncertainty bands in the x (solid line) and y (solid bold line) directions in a Montmorillonite sample with local diffusion coefficients Dil = D0 = 2.3 × 10−9 m2/s. No diffusion in the pore space between edges.

and their circularity, defined as 4πAP−2, where A is the particle area and P its perimeter. This clay sample serves as an example to demonstrate the KDPT approach capabilities in capturing the influence of pore connectivity and local diffusion properties on the average (sample scale) diffusive flux through the sample.4 In principle, for different rock types, the KDPT algorithm could easily be applied to simulate diffusion in more realistic heterogeneous domains, such as those, for example, derived from autoradiographies of real rock samples.18 As an extension to Churakov and Gimmi,4 where only four different realizations of the structure map were generated for the random walk simulations, we here propose a method to account for the stochastic nature of the smectite particles’ geometrical properties within the framework of the KDPT algorithm. The method stems from the assumption that in a compacted clay structure the diffusion of solute particles occurs mainly in the interlayers9 which, as said before, are made up of two water molecule layers. Thus, the diffusion process can be approximated as being mainly one-dimensional, at least until the tracer particle reaches the edge of the smectite particles entering the pore space between them. We further assume that the solute particle residence time within this pore space is negligible: thus, as soon as a tracer particle exits the interlayer of a smectite particle, it enters a new one in a different particle. The solute particle could in principle also diffuse back into the starting interlayer: this effect is automatically accounted for in our modeling, on average. In fact, when the new interlayer is entered, back-diffusion is allowed, although, in general, in a different direction: this minor inaccuracy occurs at different positions and orientations, thus balancing out in the ensemble of simulations. Operatively, the simulation of diffusion proceeds as follows. Initially, we consider a water tracer particle (e.g., HTO) which does not interact with the solid surfaces (ls→t = 0 and lt→s = 0)

3. RESULTS In this section, we develop the KDPT method for the process of diffusion of a solute in clay structures at the clay structure map scale.2 Due to the complexity of the structure of real clay samples and to the difficulty of obtaining direct visualizations of the pore structure at this scale by tomographic methods, real, detailed structure maps are not yet available, although some new imaging techniques, such as the focused ion beam nanotomography (FIB-nt),26 seem capable of providing promising results in the near future. Current simulation-based investigations, then, rely on artificially generated structure maps, which fulfill some average geometrical and physical properties of a real sample. In practice, the maps are built on the basis of macroscopic data, e.g., total porosity, mineralogical composition of the clay sample under analysis, average interlayer distances, particle size distributions, etc., thus retaining the main features of the pore structure of typical composite clay materials.4 The actual generation of clay structure maps accounting for the macroscopic data related to a given clay rock sample lies beyond the scope of this work. Rather, with no loss of generality, here we shall refer directly to the simplified twodimensional compacted montmorillonite structure map type provided in Churakov and Gimmi,4 which is made up of an ensemble of smectite particles of different sizes and orientations. The smectite particles, in turn, are composed of stacks of parallel clay layers, separated by interlayers whose width corresponds to two water molecules. All the relevant physical and geometrical macroscopic data of the structure map are reported in Churakov and Gimmi.4 It is worth noticing that the geometrical arrangement of the smectite particles within the structural map depends on three stochastic variables: the equivalent smectite particle radius, their angle of orientation, 18513

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

Figure 2. Average diffusion coefficients of a tracer in the x (solid line) and y (solid bold line) directions as a function of the transfer rate between the pore spaces between edges and the interlayers in a Montmorillonite sample with homogeneous local diffusion coefficients Dil = Ded = D0 = 2.3 × 10−9 m2/s.

by the tracer particle. Again, diffusion is simulated within the new one-dimensional array. Thus, each tracer particle diffuses within a different realization of the structure map. The procedure is stopped when the simulation time reaches a predefined time horizon T = 2 × 10−4 s or the particle exits the domain. Figure 1 shows the mean squared displacements ⟨x2⟩ and ⟨y2⟩ in the x (solid line) and y (solid bold line) directions computed by simulating M = 2000 tracer particles. The dotted lines represent the ±1σ uncertainty bands, where σ is the standard deviation of the sample distribution of the squared displacements of each simulation. It is worth noticing that this uncertainty accounts for the stochasticity of both the random walk, i.e., of the transition times in the KDPT scheme, and the pore structure. The estimates of the sample scale diffusion coefficients in the x and y directions, obtained by linear regression of the mean values of the squared displacements, are Dx = (⟨x2⟩/2t) = (17.88 ± 0.56) × 10−10 m2/s and Dy = (⟨y2⟩/ 2t) = (5.15 ± 0.16) × 10−10 m2/s, respectively. The uncertainty on the estimates of the diffusion coefficients is computed simply by performing a linear regression on the ±1σ uncertainty bands of the mean squared displacements. The resulting formation factors Fx = (D0/Dx) and Fy = (D0/Dy),4 which, in this case, are equal to the geometric tortuosities Gx = δsFx and Gy = δsFy because no surface effects are taken into account (i.e., the surface constrictivity is δs = 1), are Fx = 1.3 and Fy = 4.5. The sample scale diffusion coefficient D y compares satisfactorily with those obtained in Churakov and Gimmi4 for different clay structure realizations and number of random walk simulations, which assume values in the interval [4.7 × 10−10; 4.9 × 10−10] m2/s; on the other hand, the coefficient Dx is about twice as large as those found in Churakov and Gimmi,4 with values in the interval [8.7 × 10−10; 9.2 × 10−10] m2/s. To verify the validity of our model assumptions, we run a new

and with a homogeneous local diffusion coefficient in water (Dil) equal to the bulk self-diffusion coefficient of water, Dil = D0 = 2.3 × 10−9 m2/s,4 with the subscript il indicating local diffusion properties in the interlayers. The simulation of the tracer particle path in the clay structure is supposed to start at the center of a rectangular domain of size Lx × Ly = 2 × 10−6 m × 1.5 × 10−6 m. The equivalent radius Req and the orientation (with respect to a fixed coordinate system xy) of a onedimensional array representing the interlayer are then sampled from two Gaussian distributions, N(10−8 m, 4 × 10−9 m) and N(9°, 30°), respectively, inferred from data reported in Churakov and Gimmi.4 The distribution of the equivalent radius is truncated for Req < 0. The equivalent radius is then discretized in Nz = Req/Δz cells, where Δz = 8 × 10−10 m is the cell width. As verified by the authors, but not shown here for brevity’s sake, larger values of the cell size would introduce unacceptable levels of numerical dispersion in the results. Because diffusion is, in this setup, mainly one-dimensional, the third stochastic property of the smectite particles, i.e., circularity, can be here neglected. The array is centered around the initial location of the water particle, and the simulation of diffusion is then carried out according to the KDPT procedure described in the previous section. During the simulations, the current position of the particle is projected along the x and y axes, and its coordinates are tracked upon discretization of time in steps of length Δt = 2 × 10−9 s. Note that the time discretization does not affect the diffusion model (that is continuous in time) but only how the particle position is recorded. If the particle undergoes more than one transition within the same time step, then an average position is assigned to that Δt. When the tracer particle exits one of the two edges of the array, then a new equivalent radius Req and a new orientation are sampled, and the corresponding new onedimensional array is centered around the last position reached 18514

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

Figure 3. Average diffusion coefficients of a tracer in the x (solid line) and y (solid bold line) directions as a function of the transfer rate between the pore spaces between edges and the interlayers in a pyrophillite sample with heterogeneous local diffusion coefficients Dil = 2.92 × 10−9 m2/s and Ded = 1.98 × 10−9 m2/s.

consistent with the comprehensive analysis of literature data presented in Bourg et al.23 Several simulations (each with M = 2000 tracer particles) are then performed in correspondence of different values of the rate λpi, which assumes values of the form 10x, with x varying in the interval [4; 10] at regularly spaced intervals of width Δx = 0.1. Figure 2 shows the resulting sample scale diffusion coefficients Dx and Dy as a function of the rate λpi. When λpi < ≈107 s−1, once the particle exits the first interlayer, it remains in the pore space between edges for a time that is, on average, longer than the simulation time horizon: the diffusion thus occurs mainly in a random pore network with uniformly distributed orientations, and the resulting average diffusion coefficients are accordingly equal to about half of the local one (as found before). In the range 107 s−1 < ≈λpi < ≈109 s−1 the average diffusion coefficients change until they reach, for λpi > ≈109 s−1, approximately the values obtained before by neglecting diffusion in the pore space: this was to be expected because the times spent by a tracer particle in the pore space between edges tend to be very small in comparison to the typical migration transition times, thus reducing the geometrical effects of this porosity on the sample scale diffusion coefficients. Although the differences are probably due to the onedimension-based representation of the migration domain, we cannot draw any conclusive remarks on the differences in the diffusion coefficients in the x direction. 3.2. Effects of the Heterogeneity in the Local (Molecular-Scale) Diffusion Properties. The KDPT-based modeling approach described above can be used to evaluate the effects of heterogeneities in the local pore diffusion coefficients on the sample scale and average diffusion coefficients due to interactions with the surface. In Churakov and Gimmi,4 molecular dynamics simulations were performed for the diffusion of water tracers in a pyrophillite sample. Local

KDPT simulation with the same settings as the previous one but sampling the orientation of the interlayers from a uniform distribution U(0°, 360°). In this case, we obtain Dx = (11.41 ± 0.37) × 10−10 m2/s and Dy = (11.43 ± 0.34) × 10−10 m2/s, in agreement with Churakov and Gimmi4 and Kärger et al.,22 which state that in this geometrical layout one has to expect average sample scale diffusion coefficients Dx and Dy equal to half of the local ones. 3.1. Effects of the Pore Space between Edges. We now investigate the extent to which the average diffusion coefficients are affected by the tracer migration in the pore space between the edges of the smectite particles. To do so, we assume that, after exiting an interlayer, a tracer particle starts migrating in the pore space with a diffusion coefficient equal to D0. Again, under the hypothesis of high compaction, it is assumed that diffusion can be described as a one-dimensional process; thus, analogously to the modeling approach adopted for representing the interlayer structure, the pore space between smectite edges is treated as a random network of straight tubes, whose lengths are distributed according to the same Gaussian distribution used for the interlayers but with orientations distributed according to U(0°, 360°). These assumptions are in agreement with the visual representation of the structural maps given in Churakov and Gimmi.4 Moreover, it is assumed that the tracer particle re-enters an interlayer at times distributed according to an exponential distribution with rate λpi [s−1]. Within the framework of the KDPT approach, this model improvement is easily achieved by the introduction of a new particle type, i.e., a soluton in pore water between edges, which (i) can undergo the additional transition of a transformation into a soluton in the interlayer with rate λpi and (ii) migrates in a different, stochastic pore network. This approach stems from a doubleporosity representation of the diffusion domain and is similar to that adopted in Cadini et al.15 Moreover, the method is 18515

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

+

−9 Figure 4. Mean squared displacements of the Na+ ion (local diffusion coefficient DNa m2/s) in the x (solid line) and y (solid bold 0 = 1.22 × 10 line) directions in a Montmorillonite sample: (a) case with no adsorption and (b) case with transition rates ls→t = 1.32 × 108 s−1 and lt→s = 108 s−1.

simpler treatment of the clay sample geometry. Now we want to show how the KDPT approach naturally allows us to include, directly at the sample scale, a description of the relaxation of a nonequilibrium system toward equilibrium due to site-specific ion sorption and/or incorporation into the solid phase. As an example, let us consider the Na+ cation, whose + diffusion coefficient in bulk water is approximately DNa 0 = 1.22 × 10−9 m2/s.24 We assume that the weak electrostatic interaction of Na+ with the smectite surfaces can be macroscopically described by a first-order, linear kinetics of the kind of eq 1.b. In this framework, the trappon population T(z,t) representing the sorbed Na+ cations evolves in time toward the equilibrium with the soluton population S(z,t) in every cell z, so that their ratio tends to be constant, i.e., (T(z,t)/ S(z,t)) → ((ls→t)/(lt→s)), as also described in Ferrara et al.20 It turns out that, when the equilibrium is reached, eqs 1.a and 1.b are equivalent to that of the advection-dispersion under the hypothesis of linear equilibrium isotherm with a retardation factor R = 1 + ((ls→t)/(lt→s)). To underlie the effects of the firstorder, linear kinetics of Na+ sorption, we make the assumption that cation diffusion occurs mainly in the interlayers, thus neglecting migration in pore water between edges (although the KDPT model would easily allow us to account for it, as shown before). Figure 4 shows the mean squared displacements ⟨x2⟩ and ⟨y2⟩ in the x (solid line) and y (solid bold line) directions computed by simulating M = 2000 particles for (a) the case with no adsorption, ls→t = 0 s−1, and (b) the case with transition rates ls→t = 1.32 × 108 s−1 and lt→s = 108 s−1. These values are chosen so that the retardation factor which results when the equilibrium is reached is equal to the theoretical one, i.e., R = 1 + ((ls→t)/(lt→s)) = 1 + ((ρbkd)/ε) = 2.32, where ρb = 1750 kg/ m3 is the bulk dry density of the clay sample;4 ε = 0.398 is the

diffusion coefficients in the interlayers and the edge regions were obtained, yielding the values Dil = 2.92 × 10−9 m2/s and Ded = 1.98 × 10−9 m2/s, where Dil > Ded due to the hydrophobic surface properties of the pyrophillite basal planes.4 The values of Dil and Ded are then assigned to the solutons in the interlayers and those in the pore space, respectively. Several simulations (each with M = 2000 tracer particles) are then performed in correspondence of the same values of the rate λpi as in the previous set of simulations. Figure 3 shows the sample scale diffusion coefficients Dx and Dy as a function of the rate λpi. Similarly to the previous set of simulations, when λpi < ≈107 s−1, the tracer particle always remains in the pore space between the edges, and the resulting average diffusion coefficients are accordingly equal to about half of the local one, i.e., Ded. In the range 107 s−1 < ≈λpi < ≈109 s−1 the average diffusion coefficients vary until they reach, for λpi > ≈109 s−1, approximately the values Dx = (22.63 ± 0.69) × 10−10 m2/s and Dy = (6.58 ± 0.21) × 10−10 m2/s (i.e., those corresponding to neglecting diffusion in the pore space). The values obtained in Churakov and Gimmi,4 i.e., Dx = 9.2 × 10−10 m2/s and Dy = 5.0 × 10−10 m2/s, lie outside the ranges of variability that we obtained, although the values of Dy are still rather similar. Note that, as expected, the anisotropy ratio Dx/Dy remains almost the same as the case of homogeneous local diffusion coefficients, 3.47 vs 3.44, because in this case it depends only on the geometry of the strucure map, which remains the same. 3.3. Effects of Chemical Interactions between the Dissolved Ions and the Solid Matrix Surfaces. So far, we have followed a general up-scaling approach similar to that of Churakov and Gimmi.4 In fact, we have plugged local diffusion coefficients carrying the information on the surface interaction processes (found, for example, by molecular dynamics simulations) into a larger-scale model which has allowed a 18516

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

Figure 5. Mean squared displacements of the Na+ ion in the x direction in a Montmorillonite sample: (a) case with no adsorption, (b) ls→t = 1.32 × 108 s−1 and lt→s = 108 s−1, (c) ls→t = 1.32 × 107 s−1 and lt→s = 107 s−1, (d) ls→t = 1.32 × 106 s−1 and lt→s = 106 s−1, and (e) ls→t = 1.32 × 105 s−1 and lt→s = 105 s−1.

total porosity of the sample;4 and kd = 3 × 10−4 m3/kg is the Na partition coefficient: the kd value is chosen as representative of those derived from several diffusion (out and through) experiments, batch sorption measurements, and geochemical calculations carried out in Van Loon et al.25 with regards to the Benken deep borehole samples. It can be noted that the mean squared displacements in the x and y directions for the case with adsorption and desorption seem linear with time, thus indicating that the equilibrium is reached very soon. The resulting sample scale diffusion coefficients are D(a) x = (9.81 ± −10 m2/s for 0.32) × 10−10 m2/s and D(a) y = (2.84 ± 0.09) × 10 (b) −10 2 case (a) and Dx = (4.00 ± 0.12) × 10 m /s and D(b) y = (1.20 ± 0.04) × 10−10 m2/s for case (b). As expected, because we can assume that the adsorption/desorption process can be (b) well described by the linear isotherm, we obtain ((D(a) x )/(Dx )) (a) (b) = 2.45 ≈ ((Dy )/(Dy )) = 2.37 ≈ R = 2.32, where the differences are due to both the uncertainties in the KDPT estimates of the square displacements and those (not quantified here) in the least-squares estimates of the diffusion coefficients. To investigate the sensitivity of the KDPT model with respect to the values of the adsorption/desorption rates, other computations with ls→t = 1.32 × 107 s−1 and lt→s = 107 s−1 (case (c)), ls→t = 1.32 × 106 s−1 and lt→s = 106 s−1 (case (d)), and ls→t = 1.32 × 105 s−1 and lt→s = 105 s−1 (case (e)) are carried out, where the parameter R = 1 + ((ls→t)/(lt→s)) is maintained constant. Figure 5 shows the mean squared displacements ⟨x2⟩ in the x direction (⟨y2⟩ is omitted for clarity) in log−log scale, computed by simulating M = 2000 tracer particles for the cases (a), (b) (described above), (c), (d), and (e). It can be noted that the mean squared displacements ⟨x2⟩ in cases (b), (c), (d), and (e) (reactive diffusion) converge toward the same curve, meaning that, for longer times, the equilibrium is reached and the sample scale diffusion coefficients tend to be D(b,c,d,e) ≈ x,y ((D(a) x,y )/R); on the other hand, for shorter times, the probability that a particle undergoes a transition becomes very small, and

the effective diffusion coefficients remain unaffected by the adsorption/desorption process, thus yielding D(c,d,e) ≈ D(a) x,y x,y , unless the transition rates are so large that the equilibrium is reached almost at the start of the simulation (case (b)). Nonequilibrium kinetics, represented in Figure 5 by that portion of the curves between the two linear behaviors, occur at times which increase as the transition rates ls→t and lt→s decrease. At these times (or, equivalently, at the corresponding space scales), the definition of an effective diffusion coefficient becomes fuzzy. As shown above, at equilibrium the ratio ((ls→t)/(lt→s)) can be easily identified as R − 1. However, for describing the nonequilibrium dynamics, the knowledge of R is not sufficient, and the values of either lt→s or ls→t must be determined. For highly sorbing cations, e.g., Eu(III), measurements of the retardation factor R are difficult to obtain due to the extremely long times they require to reach the equilibrium in compacted clays, thus requiring knowledge of both lt→s and ls→t. In these cases, to physically derive the values of the exchange parameters lt→s and ls→t, one should resort to thermodynamic and kinetic parameters for the various ions under analysis, which, in general, are not available. A typical alternative strategy stems from estimating the model parameters by best-fitting the model results to available experimental outcomes: for example, Glaus et al.9 used inverse modeling to find the parameters of a numerically solved diffusion equation, and Cadini et al.15 applied a genetic algorithm scheme to identify the exchange parameters in a KDPT model of the transport of a tracer within a double-permeability medium. Future research will be focused on directly validating our simulation approach on some experimental results of literature and, consequently, identifying the relevant model parameters. To do so, a more advanced KDPT model will be implemented for representing 3D diffusion in a fully 3D structure map, with no additional conceptual efforts. The spatial scale of the KDPT 18517

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

prediction of the equilibrium conditions for highly sorbing cations, e.g. Eu(III), which are difficult to observe experimentally, due to the extremely long times they require to reach the equilibrium in compacted clays.

model suggests that nuclear magnetic resonance (NMR) techniques should be best suited to obtain experimental measurements, but larger-scale through-diffusion experiments may be considered as well, although further processes related to the rock sample nonhomogeneity may affect the results.27 A few examples of these kinds of experimental works are: (i) the measurements of self-diffusion coefficients of water tracer in clay samples by NMR in Nakashima et al.28 or by the throughdiffusion technique in the many experimental works cited in Bourg et al.29 and (ii) the measurements of diffusion coefficients of 22Na by through-diffusion in Glaus et al.9 and Van Loon et al.25

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +39 02 2399 6355. Fax: +39 02 2399 6309. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

4. CONCLUSIONS In this paper, a particle-tracking approach based on the Kolmogorov−Dmitriev theory of branching stochastic processes has been proposed for modeling the transport of contaminants in the clay mineral porous structures. The probabilistic description of the individual processes which may occur offers an intuitive and flexible way for describing the migration of the solute particles in the transport domain, even in the presence of complex geometries and/or physical/ chemical interaction with the solid-phase surfaces. The method has been applied to estimate the average diffusion coefficients of water and ions in a two-dimensional compacted clay sample at the submicrometer scale. The intrinsic geometrical properties of the compacted clay sample have suggested a decomposition of the original two-dimensional domain in a series of interconnected one-dimensional subdomains representing the single interlayers. This has allowed us to (i) simplify the structure map representation within the particle-tracking algorithm and (ii) automatically account for the inherent uncertainty of the structure map. The estimated diffusion coefficients of a water tracer in both the compacted montmorillonite and pyrophillite samples compare satisfactorily with those estimated by a different particletracking method in Churakov and Gimmi,4 although the component in the x-direction is somewhat larger than expected. No conclusive explanations have been found for this issue, although we suspect that the differences are probably due to the strictly two-dimension-based representation of the migration domain. However, it is worth mentioning that the main purpose of this work is that of demonstrating the applicability and the flexibility of the KDPT approach and that the extension to a full two-, or even three-, dimensional domain would be conceptually straightforward. The KDPT scheme has also allowed to easily include a double-porosity representation of the particle migration, to account for diffusion within the pore space between the edges of the clay particles. The average diffusion coefficients, estimated as a function of the exchange rate between the pore space between the edges and the interlayers, have been shown to assume values consistent with the physical intuition. Finally, the KDPT has been demonstrated to be capable of easily modeling chemical interactions of an ion with the clay solid matrix surfaces in the form of first-order, linear kinetics with constant adsorption and desorption rates. The methodology has been applied to an illustrative example of diffusion of the weakly sorbing ion Na+ within the same clay sample used as a reference in the previous case studies, showing its capability of capturing the nonequilibrium transients as well as the equilibrium behavior. The satisfactory results provided by the KDPT suggest, among others, its potential use for the

REFERENCES

(1) Chapman, N. A.; McKinley, I. G. The Geological Disposal of Nuclear Waste; Wiley: U.K., 1987. (2) Altmann, S.; Tournassat, C.; Goutelard, F.; Parneix, J.-C.; Gimmi, T.; Maes, N. Diffusion-Driven Transport in Clayrock Formations. Appl. Geochem. 2012, 27, 463−478. (3) Gimmi, T.; Kosakowski, G. How Mobile are Sorbed Cations in Clays and Clay Rocks? Environ. Sci. Technol. 2011, 45, 1443−1449. (4) Churakov, S. V.; Gimmi, T. Up-Scaling of Molecular Diffusion Coefficients in Clays: a Two-Step Approach. J. Phys. Chem. C 2011, 115, 6703−6714. (5) ONDRAF/NIRAS. Safety Assessment and Feasibility Interim Report 2 (Safir 2), ONDRAF/NIRAS 2001-05, 2001 (6) Nagra. Project Opalinus Clay, Demonstration of Disposal Feasibility for Spent Fuel, Vitrified High-Level Waste and Long-Lived Intermediatelevel Waste, Safety Report. Technical Report 02-05, 2002. (7) Andra. Dossier 2005: CDROM Argile − Plaquettes, Synthese, Tomes, Referentiels et Glossaire. www.andra.fr (reference #298), 2005. (8) Robinet, J. C.; Coelho, D.; Altmann, S. Project CatClay Deliverable (D-N°:1−1), State of the art on cation dif f usion in clay mineral systems, 2011. (9) Glaus, M. A.; Baeyens, B.; Bradbury, M. H.; Jakob, A.; Van Loon, L. R.; Yaroshchuk, A. Diffusion of 22Na and 85Sr in Montmorillonite: Evidence of Interlayer Diffusion Being the Dominant Pathway at High Compaction. Environ. Sci. Technol. 2007, 41, 478−485. (10) Kinzelbach, W. The Random Walk Method in Pollutant Transport Simulation. In Advances in Analytical and Numerical Groundwater Flow and Quality Modelling; NATO ASI Series C; Springer: Netherlands, 1987; pp 227−246. (11) Uffink, G. J. M. Analysis of dispersion by the random walk method. Ph.D. thesis, Delft University of Technology, 1990. (12) Delay, F.; Ackerer, P.; Danquigny, C. Solution of Solute Transport in Porous or Fractured Formations by Random Walk Particle Tracking: A Review. Vadose Zone J. 2005, 4, 360−379. (13) Kolmogorov, A. N.; Dmitriev, N. A. Branching Stochastic Processes. Dokl. Akad. Nauk SSSR 1947, 56, 5−8. (14) Marseguerra, M.; Zio, E. Modelling the Transport of Contaminants in Groundwater as a Branching Stochastic Process. Ann. Nucl. Energy 1997, 24, 325−644. (15) Cadini, F.; Bertoli, I.; De Sanctis, J.; Zio, E. A Novel Particle Tracking Scheme for Modeling Contaminant Transport in a DualContinua Fractured Medium. Water Resour. Res. 2012, 48, WR011694. (16) Cadini, F.; Bertoli, I.; De Sanctis, J.; Zio, E. Monte Carlo Simulation of Radionuclide Migration in Fractured Rock for the Performance Assessment of Radioactive Waste Repositories. Reliab. Eng. Syst. Safety 2013, 111, 241−247. (17) Cadini, F.; Bertoli, I.; De Sanctis, J.; Zio, E. Upscaling of a DualPermeability Monte Carlo Simulation Model for Contaminant Transport in Fractured Networks by Genetic Algorithm Parameter Identification. Stochastic Environ. Res. Risk Assess. 2013, 27, 505−516. (18) Delay, F.; Porel, G.; Sardini, P. Modelling Diffusion in a Heterogeneous Rock Matrix with a Time-Domain Lagrangian Method and an Inversion Procedure. C. R. Geosci. 2002, 334, 967−973. (19) Painter, S.; Cvetkovic, V.; Mancillas, J.; Pensado, O. Time Domain Particle Tracking Methods for Simulating Transport with 18518

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519

The Journal of Physical Chemistry C

Article

Retention and First-Order Transformation. Water Resour. Res. 2008, 44, WR005944. (20) Ferrara, A.; Marseguerra, M.; Zio, E. A Comparison between the Advection-Dispersion and the Komogorov−Dmitriev Model for Groundwater Contaminant Transport. Ann. Nucl. Energy 1999, 26, 1083−1096. (21) Bear, J. Hydraulics of groundwater; McGraw-Hill: New York, NY, 1979. (22) Kärger, J.; Ruthven, D. M.; Theodorou, D. N. Diffusion in Nanoporous Materials; Wiley-VCH: Weinheim, Germany, 2012. (23) Bourg, I. C.; Bourg, A. C. M.; Sposito, G. Modeling Diffusion and Adsorption in Compacted Bentonite: A Critical Review. J. Contam. Hydrol. 2003, 61, 293−302. (24) Lee, S. H.; Rasaiah, J. C. Molecular-Dynamics Simulation of Ion Mobility. 2. Alkali-Metal and Halide-Ions Using the SPC/E Model for Water at 25 °C. J. Phys. Chem. 1996, 100, 1420−1425. (25) Van Loon, L. R.; Baeyens, B.; Bradbury, M. H. Diffusion and Retention of Sodium and Strontium in Opalinus Clay: Comparison of Sorption Data from Diffusion and Batch Sorption Measurements, and Geochemical Calculations. Appl. Geochem. 2005, 20, 2351−2363. (26) Holzer, L.; Münch, B. Toward Reproducible Three-Dimensional Microstructure Analysis of Granular Materials and Complex Suspensions. Microsc. Microanal. 2009, 15, 130−46. (27) Gonzalez Sanchez, F.; Gimmi, T.; Juranyi, F.; Van Loon, L.; Diamond, L. W. Linking the Diffusion of Water in Compacted Clays at Two Different Time Scales: Tracer Through-Diffusion and Quasielastic Neutron Scattering. Environ. Sci. Technol. 2009, 43, 3487−3493. (28) Nakashima, Y.; Mitsumori, F.; Nakashima, S.; Takahashi, M. Measurement of Self-Diffusion Coefficients of Water in Smectite by Stimulated Echo 1H Nuclear Magnetic Resonance Imaging. Appl. Clay Sci. 1999, 14, 59−68. (29) Bourg, I. C.; Sposito, G.; Bourg, A. C. M. Tracer Diffusion in Compacted, Water-Saturated Bentonite. Clays Clay Miner. 2006, 54, 363−374.

18519

dx.doi.org/10.1021/jp405520v | J. Phys. Chem. C 2013, 117, 18510−18519