Qingrong
Xiong
*a,
Claudia
Joseph
b,
Katja
Schmeide
c and
Andrey P.
Jivkov
a
aResearch Centre for Radwaste & Decommissioning and Modelling & Simulation Centre, Dalton Nuclear Institute, The University of Manchester, UK. E-mail: qingrong.xiong@manchester.ac.uk
bGlenn T. Seaborg Institute, Physical & Life Sciences Directorate, Lawrence Livermore National Laboratory, L-231, P.O. Box 808, Livermore, CA 94550, USA
cHelmholtz-Zentrum Dresden-Rossendorf, Institute of Resource Ecology, Bautzner Landstr. 400, 01328 Dresden, Germany
First published on 26th October 2015
Compacted clays are considered as excellent candidates for barriers to radionuclide transport in future repositories for nuclear waste due to their very low hydraulic permeability. Diffusion is the dominant transport mechanism, controlled by a nano-scale pore system. Assessment of the clays' long-term containment function requires adequate modelling of such pore systems and their evolution. Existing characterisation techniques do not provide complete pore space information for effective modelling, such as pore and throat size distributions and connectivity. Special network models for reactive transport are proposed here using the complimentary character of the pore space and the solid phase. This balances the insufficient characterisation information and provides the means for future mechanical–physical–chemical coupling. The anisotropy and heterogeneity of clays is represented using different length parameters and percentage of pores in different directions. Resulting networks are described as mathematical graphs with efficient discrete calculus formulation of transport. Opalinus Clay (OPA) is chosen as an example. Experimental data for the tritiated water (HTO) and U(VI) diffusion through OPA are presented. Calculated diffusion coefficients of HTO and uranium species are within the ranges of the experimentally determined data in different clay directions. This verifies the proposed pore network model and validates that uranium complexes are diffusing as neutral species in OPA. In the case of U(VI) diffusion the method is extended to account for sorption and convection. Rather than changing pore radii by coarse grained mathematical formula, physical sorption is simulated in each pore, which is more accurate and realistic.
Pore network models (PNM) are amongst the appealing approaches that provide a suitable description of mutable pore space structures. Such models have been used to describe conservative as well as reactive transport in saturated7 or multiphase and unsaturated porous media.8,9 Moreover, their use has been expanded to complement the application of percolation theory to porous media10 and to simulate the effect of microbial growth on substrate transport and system permeability.3 The pore space is approximated by a set of sites and a set of bonds connecting some of the sites. The sites positions may form a regular lattice (regular PNM) or reflect known length variations (irregular PNM). Mass transport is allowed through the system of bonds.
Conceptually, PNM are scale indifferent, i.e. they can be applied to any length-scale interval where the structure of the pore space has been experimentally characterised. For example, if a particular experimental technique allows for characterising pore features of sizes between 0.1 μm and 10 μm, then the corresponding PNM will take these into account. The elements of PNM, sites and bonds, are abstract and can be related to the measurable features in different ways depending on the available information.
Historically, pores were related to the sites and pore throats were related to the bonds of regular PNM.11–13 If such correspondence is to be statistically representative of the material modelled, then sufficiently rich experimental information is required: shape and size distribution of pores and throats, as well as the pore coordination spectrum, i.e. percentages of pores coordinated by different numbers of throats.13,14 These can be obtained from structures with distinguishable pores and pore throats. One established way is analysing 3D images from synchrotron X-ray tomography.15,16
It is worth mentioning, that 3D images can be used directly to construct irregular PNM.17 These allow for validating physical assumptions for flow simulations using 4D imaging of mass transport. However, such irregular PNM are sample-specific and potentially not statistically representative. A regular PNM constructed within a larger volume allows for capturing statistics from a number of imaged samples, i.e. improved statistical correspondence, and for calculating transport at distances closer to the engineering application length scale.
X-ray computed tomography is presently limited to resolving macro-porosity (features larger than 50 nm). However, the resolution of this experimental technique is not sufficient to identify throats for systems dominated by meso-porosity (features between 5 nm and 50 nm) and micro-porosity (features smaller than 5 nm). For pore network construction, this means that pore connectivity data cannot be extracted as throats cannot be identified.
One possibility is to develop a regular PNM with pores from experimental distribution located at sites and notional throats between each pair of neighbouring pores1 based on a topologically proper lattice proposed by Jivkov et al.13 Such a construction allows for the calculation of the lattice length scale from known porosity and pore size distribution, but the pore connectivity is constant. The diffusivity of a notional throat depends on the sizes of the connected pores, the size of the solute molecules, and the interaction of the solute with the pore walls. The approach provided insights into the effects of pore size distribution, solute size, and sorption on diffusivity.1 It was noted, however, that the spatial randomness of pore sizes yielded variation of calculated diffusion coefficient larger than experimentally measured variability. This indicated that in the real material there was a strong connectivity effect, missed by the model due to lack of detailed information. In addition, diffusivity anisotropy emerging from material texture (preferred mineral orientations) could not be captured.
Another possibility is to assume variable connectivity.18 While this was more realistic in terms of pore space topology, the observable data were insufficient to establish rigorously a lattice length scale. With variable length scales, improved agreement with the variability in experimentally measured diffusivity was demonstrated, including diffusion anisotropy in clay.
The two approaches to tackle incomplete pore space information, namely predefined connectivity for calculable length scale,1 and undetermined length scale for realistic connectivity,18 suffer from the lack of an additional constraint. This cannot be found within the pore space information and requires the consideration of the solid phase structure, e.g., the shape and size distribution of mineral grains. A methodology for incorporating the structure of the solid phase is proposed here. It improves substantially the realism of the constructed PNM, both in terms of geometry and topology. The proposal has added benefit that the constructed PNM can be paired directly to lattice models of the solid phase developed for analysis of damage evolution via micro-cracking.19–21 This facilitates mechanistic investigations of combined mechanical–thermal–chemical–biological effects on diffusivity, which will be discussed further in the paper.
Typical pore network models idealise the pore space with simple geometries, i.e. throats are assumed to have square or circular cross sections. This makes it difficult to simulate processes like sorption, which may change pore geometries by physical clogging. Recent advances allowed to model pore geometries with irregular cross sectional shapes.8 The influence of cross sectional shapes on diffusion and sorption will be analysed in this work.
For safety assessment of repositories for high-level nuclear waste, in particular the migration behaviour of radioisotopes with long half-lives (e.g., 238U, τ1/2 = 4.468 × 109 a; main fraction of spent nuclear fuel) is in the focus of research. Opalinus Clay (OPA) is discussed as potential host rock for nuclear waste repositories, where it is considered as one of the barriers to radionuclide transport, because it exhibits very low hydraulic conductivity.22 This makes molecular diffusion the main transport mechanism of radionuclides through the barrier. In the present work, OPA is selected to demonstrate the development of the new pore network model. The pore space information for model construction is adopted from Keller et al.23,24 and Nagra.22 Validation data on macroscopic diffusivities come from recent literatures.25–28 Although OPA has a very low permeability, the effects of velocities on mass transport are also analysed to demonstrate how the methodology will be applied to other porous media. In the diffusion experiments so far the neutral Ca2UO2(CO3)3(aq) complex was assumed to be the main diffusing species through OPA,25 because it was spectroscopically verified to be the dominant U(VI) species in the experimental source reservoir solution. The model described in the present study shall clarify whether U(VI) diffuses as cationic, anionic, or neutral species through OPA.
While the focus of this work is on the reactive transport of radionuclides through Opalinus Clay, the modelling approach proposed is applicable to analysis of transport phenomena in other porous media, where the pore space cannot be fully characterised by current experimental techniques. Examples include reactive transport of species in micro-porous and meso-porous materials which may not be easily studied by classical in-diffusion or through-diffusion methods, biomass growth in bioengineering materials and dissolution–precipitation processes in carbon capture and storage.
The material is reduced to a mathematical graph by placing particles or grains at cell centres (interiors). This is illustrated in Fig. 1(a) for cells with equal distances between the three pairs of square faces, a setup used in the previous works. In order to represent texture anisotropy, three different length parameters, S1 [m], S2 [m], and S3 [m], are used to measure the distances between the square faces in directions (1, 0, 0), (0, 1, 0), and (0, 0, 1), respectively. All other distances are then dependent on the three length parameters, e.g., the volume of each cell is Vc = S1S2S3/2 [m3].
In an assembly of Nc [—] cells, each cell is assigned to a particle with different radius, ri [m]. The assignment of particles terminates when their accumulated volume fraction equals to the experimentally measured value. From this requirement the volume of a cell assembly is calculated by:
(1) |
Each particle (site, cell centre) is connected to its 14 neighbouring particles by bonds representing relative deformations between cells, illustrated in Fig. 1(b). This graph can be analysed very efficiently using discrete exterior calculus to investigate the response of the material to loading, including micro-crack generation, competition, and coalescence.32 The presence of pores in cell neighbourhoods is reflected in pore-size-dependent reduction of bond strengths. Further, the assigned particle sizes may also contribute to spatial variation of bond properties within the model. However, for the construction of the PNM the particular particle sizes do not play a role; only the length parameter calculated by eqn (1) is of essence.
The pore network construction is as follows. Firstly, a new skeleton is formed using sites at the centres of cell faces and bonds between neighbouring faces. This is illustrated in Fig. 1(a) on a single cell with S1 = S2 = S3 for clarity. The bonds represent potential transport pathways, i.e. they show the positions where elongated pores are allowed to reside. Note that there are no bonds cutting through cells, hence the pore system can follow paths around particles. The elongated pores can be assigned with different cross sectional shapes, such as circular, equilateral, and square. Since pores reside on bonds, the sites at face centres represent pore junctions – volume-less containers redirecting mass transport.
Notably, this skeleton allows for pores with maximum coordination of 8 or 12 depending on whether the bond connects square to hexagonal or hexagonal to hexagonal boundary. Importantly, the length parameters of the skeleton, determining the lengths of possible pores, are now dictated by the solid phase structure.
(2) |
(3) |
Secondly, micro-pores are assigned to bonds not already occupied by meso-pores. The micro-pores are also assumed to be of different shape and each new pore assignment increments the pore volume according to selected cross section area and bond length. The process terminates when the total pore volume fraction, from meso- and micro-pores, attains the experimental porosity. Notably, micro-pores are allowed at bonds in all directions, including bonds in the bedding direction not occupied by meso-pores, and more than one micro-pore can be assigned to a bond if required to reach prescribed porosity, illustration in Fig. 2(c). A set of micro-pores assigned to a bond can be considered as providing an “effective” transport between two pore junctions, which is necessary to create transport paths in the otherwise non-percolating meso-pore system. Importantly, they act as sieves for transport of species with different molecular sizes.
An example of a PNM after assigning meso-pores and micro-pores with different pore shape sets is shown in Fig. 2(b). The micro-pore sets are presented as bundles with scaled-up cross-sections and the network is constructed with equal length parameters in the bedding and out-of bedding directions for improved illustration. In the simulations different length parameters in bedding and out-of bedding directions were used to represent texture.
Regarding the pore space, we make use of the 3D analysis of OPA reported by Keller et al.23,24 and Nagra.22 The pore space was divided into meso-pores and micro-pores based on the resolution of different experimental techniques in our work. The meso-pores were resolved by Focused Ion Beam nano-tomography (FIB-nt) with elongated pores in the bedding plane. These pores had approximately cylindrical shapes with diameters >10 nm and were largely isolated (i.e. solute species cannot percolate through the sample). The meso-pores occupied approximately 1.8 vol% of the sample, i.e. their volume density (porosity) was θmes = 0.018. The micro-pores were resolved by N2 adsorption with pore sizes <10 nm, but were unresolved by FIB-nt. The analysis of N2 adsorption gave a total physical porosity of around 11.5 vol%, of which 9.7 vol% were micro-pores. The volume density of micro-pores is thus θmic = 0.097.23 However, micro-pores narrower than 2 nm cannot be probed by N2 adsorption. For OPA from Mont Terri, most estimates of the total porosity are on the order of 16 vol%.25,36,37 In this work, we assume that micro-pores smaller than 2 nm occupy 4.5 vol% and follow the same distribution as pores occupying 9.7 vol%. These measurements show that the fine-grained clay matrix contains an extensive micro-pore network. The meso-pore system can be seen as part of the effective “flow relevant” porosity that surrounds the clay mineral particles and is interconnected by micro-pores.
The two measurements were combined into a single “cumulative pore volume fraction – pore radius” curve given in Fig. 3(a).18,23 This shows that the geometrically resolved meso-pores (symbols) and the adsorption-based curve coincide in their common length interval, with meso-porosity θmes = 0.018 and total porosity θ = 0.16. The radii are derived with the assumption of cylindrical pore shapes truthful to experimental evidence for meso-pores23 and approximate for micro-pores. For the purpose of model construction the experimental distribution from Fig. 3(a) is re-evaluated as cumulative probability separately for meso- and micro-pores. These are shown in Fig. 3(b) and (c), respectively. The use of the cumulative probability is to produce sizes from the measured distribution in Fig. 3(a) using a generator for uniformly distributed numbers 0 ≤ p < 1. This is illustrated by the red arrows on Fig. 3(b) and (c); more details on statistics are given by Jivkov and Xiong.18 For pores with other shapes such as square, the radii are converted to length of side with equal pore volume.
Fig. 3 Microstructure characteristics of OPA: (a) cumulative pore volume fraction versus pore size determined by FIB-nt and N2 adsorption analyses;23 (b) cumulative distribution of meso-pore sizes; (c) cumulative distribution of micro-pore sizes with assumed probability density from ref. 18; (d) cumulative distribution of carbonate grain sizes.24 |
Regarding the structure of the solid phase, Keller et al.24 reported 18 vol% of non-porous carbonates and 17 vol% of non-porous quartz in the sample. The sizes of carbonates ranged from 100 to 300 nm, while the size distribution of quartz was not measured. The reported data of carbonate were converted to the cumulative probability of its grain sizes shown in Fig. 3(d) for model construction. The same cumulative probability is assumed for the non-porous quartz, due to the lack of more accurate data, and used to generate particle sizes of both types in a way similar to pore size generation.
(4) |
To verify the pore network model, two diffusion species are considered in this work. One is HTO, a neutral solute with molecular size of rA = 0.1 nm.18 The other is a neutral U(VI) complex, Ca2UO2(CO3)3, which dominates the U(VI) speciation in OPA pore water at room temperature to nearly 100%.40 The molecular size of this complex is rA = 0.524 nm.41 These two species are selected for comparison with experimental data on their diffusivity in OPA. The free molecular diffusion coefficients calculated with the Einstein–Stokes equation at room temperature, are D0 = 2.24 × 10−9 m2 s−1 for HTO42 and D0 = 4.66 × 10−10 m2 s−1 for U(VI),43 respectively.
(5) |
(6) |
(7) |
(8) |
(9) |
(10) |
(11) |
A temporal and spatial discretization of above equations is used to describe mass transport through a pore
(12) |
(13) |
Linear equilibrium adsorption is described by:
sij = ka,ijcij, | (14) |
Then eqn (13) becomes
(15) |
For nonlinear equilibrium adsorption the Freundlich isotherm48 is used:
sij = acbij, | (16) |
(17) |
The incidence matrix describes both the topological structure of the system, which can be used for connectivity analysis, and the derivative of a discrete function defined on the nodes.49 Specifically, a concentration field on the graph nodes, vector of dimension N, has a gradient given by the matrix product , which is a discrete field on the edges, vector of dimension E. As importantly, the transpose of the incidence matrix describes the derivative of a discrete function defined on the edges. Specifically, a mass flow field on the graph edges, a vector of dimension E, has a gradient given by , which is a discrete field on the nodes, a vector of dimension N.
The components of the discrete field of edge flows, , are calculated by eqn (15) or (17), i.e. each component is the edge gradient scaled by an edge weight given by We = (qij + DijAij/Lij)/(1 + kaβ) or We = (qij + DijAij/Lij)/(1 + acb−1ijβ) for linear and nonlinear sorption at pore scale, respectively. Clearly, the edge weight depends on the assigned pore radius, velocity in the pore, selected length scale, radius of diffusing species, and sorption coefficient. For convenience the edge weights are arranged in a diagonal matrix, , of dimensions E × E, where the element in row and column e is the weight of edge e. Thus, the discrete flow field is given by .
The mass transport through the graph is governed by mass conservation at nodes, i.e. the rate of change of diffusing species mass at the nodes equals the resultant mass flows through the nodes. The latter is the algebraic sum of the flows through adjacent edges, which is the divergence of the flow field given by . Thus, the process of mass transport through a weighted graph is formulated mathematically by
(18) |
The formulation of flow via discrete analysis on graphs is highly beneficial for incorporating system evolution compared to traditional finite element or finite difference formulations. Conceptually, it shows the intrinsic relation between the pore space topological structure and transport, both described via the incidence matrix. Practically, it differentiates clearly between effects of connectivity, captured by the incidence matrix , and geometry and physics, captured by the edge weight matrix . This facilitates the incorporation of potential pore space changing mechanisms in a computationally effective way. For example, local pore geometrical changes arising from corrosion, sorption, bacterial film growth, deformation, etc., will affect corresponding coefficients of only but leave constant. Material failures, such as micro-cracks, during evolution present a more complex case. These change the connectivity, i.e. they alter , and introduce geometry and physics for the new edges, i.e. modify .
Immutable pore systems have been analysed for the calculation of macroscopic diffusivity. Comparing the calculated results with experimental diffusion coefficients reported in the literature, the realism of the model construction in terms of initial pore space geometry and topology will be verified. In this work, the effects of pore shapes, sorption coefficient, and velocities are analysed based on previous work.
Fig. 4 Schematic depiction of the experimental set-up used in the diffusion experiments (based on Van Loon et al.50). |
At first, HTO through- and out-diffusion experiments were performed according to the procedure described by Van Loon et al.50 in order to determine values for the total porosity (θ) of the clay samples (c0(HTO) = 1000 Bq mL−1). After that, the U(VI) in-diffusion was studied (c0(233U(VI)) = 1 × 10−6 mol L−1). After 3 months, the diffusion experiments were stopped and the clay samples were removed from the cells. With the help of the abrasive peeling technique,52 U(VI) diffusion profiles were determined. The peeled layers were extracted for U(VI) content by 1 mol L−1 HNO3 (p.a., Roth, Karlsruhe, Germany). The tracer concentrations in the extracts were determined by liquid scintillation counting (mod. TriCarb 3100 TR, Perkin Elmer, Freiburg, Germany). All experimental results were evaluated using the commercial software COMSOL Multiphysics: finite-element software package 3.5a (2008) and 4.2a (2011). A detailed description of all diffusion experiments can be found in Joseph et al.25,27
Zero concentration of diffusing species is taken as initial conditions in all nodes. The boundary conditions reflect a particular experimental setup, where concentrations c0 and c1, are prescribed on two opposite boundaries, while the remaining four boundaries do not permit flux. The selection of the two boundaries depends on the macroscopic transport properties being analysed. Specifically, the boundary conditions used to calculate the macroscopic transport properties parallel to bedding direction (1, 0, 0), D1, are: prescribed concentration c0, in all nodes on plane X1 = 0; prescribed concentration c1 in all nodes on plane X1 = 20S1; zero flux through all nodes on planes X2 = 0, X2 = 20S2, X3 = 0, X3 = 20S3. For calculating the macroscopic transport properties perpendicular to the bedding direction, say (0, 1, 0), D2, the boundary conditions are: prescribed concentration c0 in all nodes on plane X2 = 0; prescribed concentration c1 in all nodes on plane X2 = 20S2; zero flux through all nodes on planes X1 = 0, X1 = 20S1, X3 = 0, X3 = 20S3.
One example of diffusion in the bedding direction for a system with unequal length parameters and circular cross sectional pores is shown in Fig. 5 by the solute concentration profiles at select times. The prescribed concentration c0 on plane X1 is 1 mol m−3 and the prescribed concentration c1 on plane X0 is 0. These are snapshots taken from the whole history obtained by numerical simulation based on eqn (4). The profiles are shown on the cellular assembly, rather than on the graph representing the pore network, for clarity.
Fig. 5 Concentration profiles of HTO [mol m−3] at time t = 0.1 s (a), t = 0.5 s (b), and t = 1 s (c). The pores are assumed with circular cross section. |
The diffusion of tracer molecules in clay has shown that mainly neutral and cationic species diffuse in the interlayer water, diffuse layer water and the free pore water.53 Experiments reported that the accessible porosity for HTO25,26,37,54 agrees very well with total porosity in OPA,22 which implies that almost all the pores in OPA are available. The total porosity has been used to simulate the effective diffusion coefficients of HTO in this work. In Fig. 6, the effective diffusivities of HTO and U(VI) parallel and perpendicular to the bedding, D1 and D2, are shown to increase as the average shape factor G increases. The calculated effective diffusivity of HTO is in the following ranges: D1 = 5.41 × 10−11 ∼ 17.10 × 10−11 m2 s−1; D2 = 1.58 × 10−11 ∼ 5.03 × 10−11 m2 s−1 (see Fig. 6(a)). Reported experimentally obtained values for HTO diffusion in OPA shown in Fig. 6(a) are in the range of D1 = (5.4 ± 0.4) × 10−11 m2 s−126 and D2 = (1.48 ± 0.03) × 10−11 m2 s−1,25 (1.26 ± 0.09) × 10−11 m2 s−1,25 (1.23 ± 0.06) × 10−11 m2 s−1,26 (1.67 ± 0.12) × 10−11 m2 s−1,27 (1.33 ± 0.11) × 10−11 m2 s−127 determined for different bore core samples. In OPA, the accessible porosity is different for the various species resulting in different effective diffusivities. The increased diffusion of cations is explained either by interlayer or by diffusion in the electrical double layer37,55 that surrounds the negatively charged clay surface and contains an excess of cations.54,56–61 The diminished diffusion of anions is due to their limited access to the interlayer space. This is caused by the overlay of double layers. As a result the electric potential becomes considerably large.62
If in the model the diffusing U(VI) species is neutral, then the same initial accessible porosity as for HTO can be assumed. The simulated effective diffusivity of U(VI) in the bedding direction, D1, (cf.Fig. 6(b)) is in the range of 8.04 × 10−12 ∼ 2.56 × 10−11 m2 s−1, and in the perpendicular direction, D2, is from 1.98 × 10−12 to 7.26 × 10−12 m2 s−1. Experimental data for U(VI) in the perpendicular direction are 1.5 × 10−12 ∼ 2.3 × 10−12 m2 s−1.25,27 The simulated results for HTO and U(VI) are in good agreement with experimental investigations when average shape factor, G, is equal to 0.025 (cf.Fig. 6). In OPA, three types of water accessible for diffusion can be distinguished. The solute, for instance a U(VI) species, enters freely the larger pores, since it resides in an electrically neutral solution, the so-called free pore-water. On the net negatively charged clay surface, the water turns into a diffuse double layer. The charge deficiency can be balanced by exchangeable cations, anions are repulsed from the surface. The third type is interlayer water in the smectite fraction of OPA. It is considered to contain exchangeable cations only, which behave chemically and physically comparable to the surface complexed cations. The types of water or pore space accessible for diffusion depends on the diffusing species.22,37,63,64 Specifically, diffusion of cations is enhanced, and that of anions diminished, relative to a neutral tracer such as HTO.37 The here discussed effective diffusivity excludes sorption effects of the U(VI) species, Ca2UO2(CO3)3(aq), on OPA during their migration, only the pure transport process is considered. Since the same approach was used for modelling HTO and U(VI) diffusion except that the radius of the diffusing species was changed, the observed difference of one order of magnitude in D2 for HTO and U(VI) can be attributed solely to size exclusion effects and not to charge effects of attracted or repulsed charged U(VI) species. It has been demonstrated in previous works,54,65 that the amount of electrical double layer water is approximately equal to half of the porosity. This means that anions cannot access almost half of the porosity, manifested in anion exclusion effect on diffusivity. Hence, the effective diffusivity of UO2(CO3)34− and U(VI) species should be different despite of their similar molecular sizes. Therefore, it can be concluded that the neutral Ca2UO2(CO3)3(aq) complex is not only the dominant U(VI) species in the OPA pore water solution in the source reservoir but also the main diffusing U(VI) species through OPA.
Overall, the model in this work can predict diffusion in OPA very well. The difference between the computational and the experimental results could be partially due to a difference between the microstructure characteristics of OPA obtained by Keller et al.23,24 used for model construction, and the clays used for experimental measurement of effective diffusion coefficients by Van Loon et al.26 and Joseph et al.25,27
Because effective diffusivities of U(VI) agree well with experimental values when the average shape factor, G, equals to 0.025, this value of G is adopted to study sorption effects of U(VI) in OPA. We assumed that linear sorption occurs at pore-scale according to experimental investigation.35 Apparent diffusivities of U(VI) species are shown in Fig. 7(b) as a function of distribution coefficients, ka. It is clear that the increase of ka yields reductions of the apparent diffusivity, D1a and D2a. Because the adsorbed mass is dependent on the value of ka, more species will be adsorbed when ka increases. As a consequence, migration of the U(VI) species retards because of their stronger interaction with the OPA pore walls. In addition, since more species are adsorbed onto the pore-wall, the cross section of pores allowed for diffusion becomes smaller and thus the accessible pore space will be reduced. Notably, the decrease of the apparent diffusivity becomes slowly as the adsorption coefficient increases. This is an illustration of the non-linear effect of the sorption in the system.
The calculated results of apparent diffusivities are in the following ranges: D1a = 8.35 × 10−14 ∼ 2.05 × 10−13 m2 s−1; D2a = 1.98 × 10−14 ∼ 4.87 × 10−14 m2 s−1, see Fig. 7(b), when ka decreases from 250 nm to 100 nm. Reported experimentally obtained values for U(VI) complexes diffusion in OPA are D1a = 9.5 × 10−14 m2 s−1,28D2a = 2.8 × 10−14 ∼ 3.4 × 10−14 m2 s−1,25,27 2.8 × 10−14 m2 s−1.28 These ranges are given in Fig. 7(b). Direct comparison shows that the simulation and the experimental results are in reasonable agreement when ka is about 200 nm. The experimental distribution coefficients kd measured by experiments were parallel to the bedding: 0.0142 m3 kg−1,28 perpendicular to the bedding: 0.025 ± 0.003 m3 kg−1,25 0.038 ± 0.005 m3 kg−1,27 0.0138 m3 kg−1.28 Since ka is the distribution coefficient, kd, normalised by surface area, and the surface area obtained by ethylene glycol monoethyl ether for OPA varies from 80 to 107 m2 g−1 (the average value is 94 m2 g−1),22 the experimental ka values range from 128 nm to 350 nm. This is very close to the simulated value of ka. The sorption of U(VI) species is therefore linear.
The pore-scale Freundlich isotherm sorption, s = acb (cf.eqn (16)), is also studied in this work. Since this sorption equation is dependent on the pore-local concentration variation of solutes, the concentration applied at the boundaries must be clarified. This is different for the analysis of effective diffusion and linear sorption where the results are independent of the prescribed concentration difference on the boundaries. We apply C0 = 0.05 mol m−3 at X0 = 0, C1 = 0 at X1 = 20S1, and no flux at other boundaries. The effects of parameters a and b on the apparent diffusivities are shown in Fig. 8. In Fig. 8(a), the increase of a yields a reduction of the apparent diffusivities of U(VI) in both directions when b = 1.5. The changes are similar with the situation of linear sorption, s = kac. In Fig. 8(b), the apparent diffusivities of U(VI) in both directions increase as b increase. The increase of apparent diffusivities becomes slowly with increasing b. This is due to the fact that the adsorbed mass becomes less as the parameter of b increases. The maximum values of apparent diffusivities will be very close to their effective diffusivities when b increases.
Fig. 9 Time history of uranium concentration at the outlet surface for different velocities, the distribution coefficient ka is 200 nm, the pore shape factor is 0.025. |
The effect of clay's texture on the transport properties has been investigated by analysing the transport of HTO and U(VI) through OPA. The anisotropy and heterogeneity of the clay can be simulated via different length parameters and percentage of pores in different directions. Different pore shapes are also taken into account in this model, which is very important for multi-phase transport. Calculated diffusion coefficients are within the ranges of experimental data. The results suggest that U(VI) diffuses as neutral species, Ca2UO2(CO3)3(aq), through OPA. Because only species in the source reservoir solution are known in experiments, the methodology proposed in this work could clarify the kind of U(VI) species (anion, cation or neutral species) diffusing in OPA. In addition, the idea has been extended to account for sorption and convection. Physical sorption is simulated in each pore, which is more accurate and realistic compared to previous coarse grained mathematical formula by changing pore radii.1
Although this work has used a particular material to illustrate the methodology and its application, the model is suitable for modelling and simulation of other micro- and meso-porous materials. The methodology proposed in this work can be used to study sorption of species that may not be easily studied by classical in-diffusion or through-diffusion methods, for example, when the species is adsorbed onto laboratory materials (e.g., cells, vessels, and tubing).
The proposed model can be potentially used for coupling with deformation and fracture of porous media, where the pore system and the solid phase, represented as dual graphs, will interact effectively. Models mimicking the coupling are possible and are subject of ongoing work. For example, our immediate plans are to study the effect of micro-cracking due to thermo-mechanical loads on the diffusion with application to long-term behaviour of OPA in a nuclear wastes repository. The goal with the concept development is to have a set of discrete models, i.e. PNM accounting for transport mechanism and lattice models of the solid phase developed for analysis of damage evolution via micro-cracking,19 that evolve concurrently at their own rates and inform each other about changes in their structures. Such a coupled problem could be of great interest in many disciplines, e.g., shale gas extraction, dams, and underground cavities in hydropower engineering. In these applications, the evolution of permeability rather than diffusivity of the medium is required, but the methodology developed in this work would be the same.
This journal is © the Owner Societies 2015 |