Study on Tritium and Iodine Species Transport through Porous Granite: A Non-Sorption Effect by Anion Exclusion

The safety of deep geological repositories is important in the disposal of high-level radioactive waste (HLW). In this study, advection–dispersion experiments were designed to build a transport model through a calibration/validation process, and the transport behavior of tritiated water (HTO) and various iodine species (iodide: I− and iodate: IO3−) was studied on a dynamic compacted granite column. Breakthrough curves (BTCs) were plotted under various flow rates (1–5 mL/min). BTCs showed that the non-sorption effect by anion exclusion was observed only in I− transport because the retardation factor (R) of I− was lower than that of HTO (R = 1). Moreover, equilibrium and nonequilibrium transport models were used and compared to identify the mobile/immobile zones in the compacted granite column. The anion exclusion effect was influenced by the immobile zones in the column. The non-sorption effect by anion exclusion (R < 1) was only observed for I− at 5.0 ± 0.2 mL/min flow rate, and a relatively higher Coulomb’s repulsive force may be caused by the smaller hydration radius of I−(3.31 Å) than that of IO3−(3.74 Å).


Introduction
It is necessary to carry out safety evaluations for deep geological disposal site of high-level radioactive waste (HLW) to evaluate the migration behavior of nuclides in surrounding rock [1]. In general, a common cause of radionuclide migration is groundwater intrusion into HLW repositories due to the failure of engineering barriers. Granite is a potential crystalline host rock in HLW geological disposal [2][3][4][5].
Considerable research efforts worldwide have increased our knowledge and understanding of how underground disposal systems function over long periods of time. Recently, several studies have investigated the advection-dispersion behavior of cationic radionuclides, such as 238 Pu, 237 Np, 137 Cs, and 90 Sr [6][7][8][9][10], to evaluate the safety of the migration of various radionuclides in deep geological environments. At present, cationic radionuclides in liquid show high adsorption ability on mineral surfaces because of the permanently negative charges by isomorphous substitution in the crystal lattice of a mineral. Migration of cationic radionuclides in host rocks would be retarded, but anionic radionuclides (i.e., 36 Cl, 99 Tc, 129 I) show exclusive or repulsive behavior by Coulomb's force [11][12][13]. Few studies on safety assessments (SA) for HLW disposal focused on the anion exclusion effect of anionic nuclide transport in granite and other host rocks [14][15][16].
In addition to transuranic radionuclides (TRUs) in HLW, the long-lived radionuclides that exhibit a relatively high mobility under deep geochemical conditions must be considered when building an HLW repository. As a primary fission product, 129 I is characterized by a long half-life (t 1/2 = 1.57 × 10 7 years), high fission yield, easy volatilization, easy migration, high radioactive toxicity, high bioavailability, and weakly adsorbed radionuclide anions [17][18][19][20]. Actually, iodine species that have various valence states (iodide (I − ), iodate (IO 3 − ), I 2 , and organic iodine) depend on pH and redox conditions in aqueous solution. Thus, 129 I is the main source of the potential risk posed not only by HLW repositories but also by nuclear accidents, such as the 2011 Fukushima Daiichi nuclear disaster. This threat was recently made obvious when 129 I was observed in the Pacific Ocean around Japan, where fast dispersion by seawater movement and deposition caused it to migrate into the marine sediments [21].
The in-house dynamic column technique and numerical analyses have been simultaneously applied to simulate the transport of different contaminated solutes under natural groundwater conditions to provide a safety assessment for environmental impact on industrial or radioactive waste storage and disposal. According to different experimental conditions, the tests are designed to the filling mode of geotechnical medium (undisturbed or broken), solute injection mode (pulse or continuous injection), and water content maintenance mode (unsaturated or aquifer zone) [22][23][24]. To understand the anion exclusion effect, a previous study injected tritiated water (HTO) and anions into the same column and expressed the "accelerated" migration of anions by comparing the experimental breakthrough curves (BTCs) with the effluent concentration [25]. By comparing the BTCs of HDO and Cl − in the dual-porosity structure of peat, McCarter obtained the anion exclusion of Cl − and found that anion exclusion weakens with increasing Cl − diffusion in pores [16]. Rao and Jain [25] determined the "accelerated" transport of [CO(CN) 6 ] − under the anion exclusion effect by analyzing the migration concentration distribution curves of HTO and [CO(CN) 6 ] 3 − in clay and established a numerical model. In addition, some studies proposed "accelerated" migration by comparing the peak position of BTCs with pore volume (PV) (if PV < 1, then accelerated transport occurs). Shukla studied the transport behavior of Cl − in soil columns of different particle sizes and indicated the anion exclusion by comparing the peak position. However, the homogeneity of the filling soils must be guaranteed to prevent the preferential flow (PV < 1) in the column.
In 1972, a two-region numerical model including the immobile/mobile zone in saturated porous soil based on diffusion double layer (DDL) theory was first developed by Krupp for the solute's pore velocity distribution formed by anion exclusion [26,27]. In addition, a physical nonequilibrium ion exclusion model was proposed and modified from DDL, and the porewater in porous soil was separated and expressed as flowing and nonflowing areas. Moreover, field experimental results by Gvirtzman [28] explained a dispersion model without ion exchange and an exchange model with ion exchange and obtained a good fitting.
In a previous work [29,30], individual and coexisting I − and IO 3 − were successfully analyzed by applying ion chromatography (IC) coupled with a convenient and rapid detection tool for inductively coupled plasma optical emission spectrometry (ICP-OES), and it has become a quick and effective method in analyzing anion speciation in environmental samples. In the present study, the advection-dispersion behavior of I − and IO 3 − in granite was studied and used to simulate their transport behavior in a micropore system. The BTCs of HTO, I − , and IO 3 − at different flow rates were plotted and compared, and the key transport parameters were estimated from STudio of ANalytical MODels (STANMOD) and HYDRUS-1D under various conditions. The R, water content in the immobile regions (θ im ), first-order mass transfer coefficient (α), and other parameters were obtained to quantify the anion exclusion effect, and the influence of flow rate and iodine speciation on anion exclusion was explored. The results will support safety assessments of the impact of "accelerated migration" of anionic radionuclides for HLW repositories.

Equilibrium Transport
The majority of studies focused on radionuclide transport in saturated porous media, such as compacted crushed granite, and the one-dimensional governing solute transport equation is: where C is radionuclide concentration in solution [Bq/mL 3 ], v is the average pore-water velocity [cm/min], t is time [min], and x is the distance [cm]. D is the dispersion coefficient of the granite media [cm 2 /min], and it can be further defined and written as where D L [cm] is the dispersivity and D * the diffusion coefficient in free water [cm 2 /min].
For the case where the adsorption equilibrium can be reached instantaneously or it takes a long time to reach the adsorption equilibrium, the distribution coefficient K d is expressed as where K d is the distribution coefficient [mL 3 /g]. In this study, we only take account of instant adsorption of radionuclides on the pores of crushed granite, and Equation (2) can be re-written as ∂C ∂t where θ is the porosity or volumetric water content [-], and ρ b is the bulk density of porous medium [mL 3 /g]. It can be combined by a retardation factor given by R = 1 + ρ b K d θ , and Equation (3) can be rewritten as In addition, the boundary conditions and initial conditions are where C 0 is inlet radionuclide concentration [Bq/mL 3 ], and L is the length of the column [cm].

Two-Region Nonequilibrium Transport
A two-region transport model was proposed and assumed that the liquid phase can be partitioned into mobile (flowing) and immobile (stagnant) regions (which are produced by anion exclusion). The solute exchange between the two liquid regions is simulated as a first-order kinetic process. The two-region solute transport model is given by where the subscripts m and im refer to the mobile and immobile liquid regions, respectively, J w = v × θ = v m × θ m is the volumetric water flux density, f represents the fraction of adsorption sites that equilibrates with the mobile liquid phase, and a is the first-order mass transfer coefficient governing the rate of solute exchange between the mobile and immobile liquid regions. θ is equal to θ m + θ im . µ l,m and µ l,im are first-order decay coefficients for the mobile and immobile liquid phases, respectively; µ s,m and µ s,im are first-order decay coefficients for the mobile and immobile adsorbed phases, respectively; r l,m and r l,im are zero-order production for the mobile and immobile liquid phases, respectively; and r s,m and r s,im are zero-order production terms for the mobile and immobile adsorbed phases, respectively.  Table 1 lists the parameters of the filled granite, the composition of the synthetic groundwater, and the size of the column. The inner diameter of the experimental column was 1.6 cm, and the length was 30 cm. A total of 85.40 ± 0.50 g crushed granite was filled, with a density of 1.45 ± 0.05 g/cm 3 and a porosity of 0.47 ± 0.03. In order to provide a stable experimental flow rate, the device was also equipped with a pressure sensor and switch drain valves.

Rocks and Liquids
A discrete granitic island surrounding Fujian Province and local synthetic GW were sampled and investigated in previous studies [6,10]. All granite samples were crushed with a grinder to a particle size of 1 mm. Before filling the soil column to crush granite samples, all granite samples were washed three times with deionized water, dried at 100 ± 10 • C for 12-48 h, and stored in a vacuum dryer for future use.

Mineral Composition and Elemental Analysis
A mineral analysis of granite in rock was conducted using X-ray diffraction (XRD) as previously described. XRD spectra and micro-polar microscopy images showed that the main minerals included quartz, plagioclase, K-feldspars, amphibole, and biotite [31,32]. In addition, the crushed granite samples (~2 g) were analyzed under scanning electron microscopy coupled with energy dispersive spectrometry (SEM-EDS) with an accelerating voltage of 20 kV and a current of 10 µA. EDS was used to analyze the major elements of the granite samples.

Pre-Equilibrium: Water Saturation in Column
In order to simulate the migration behavior of nuclides under aquifer conditions, it was necessary to saturate the crushed granite column. The valve was opened to make the groundwater flow into the column slowly at a flow rate of 3 mL/min. During the water filling, 10 mL of effluent from each column was extracted every 30 min. The concentrations of K, Ca, Mg and Na in the effluent were measured by ICP-OES. When the concentration changed within 5% of the corresponding liquid phase concentration, it indicated that water saturation had been reached.

Nonreactive Tests of HTO
In order to characterize the internal structure of the crushed granite column, the dispersion and effective porosity were obtained by convection dispersion experiments with nonreactive tracer (HTO). By regularly collecting the effluent at different flow rates, the HTO content in the effluent was analyzed and measured, and the data were used to draw the penetration curve BTC. The water circulation was divided into two stages. The first stage was the manual operation stage, and its basic process was as follows. As shown in Figure 1, set water circulation under different flow rate conditions, 3 PV (about 2100 Bq) each time, including up/down flooding. The liquid sample was collected by an automatic sampler, and a 5 mL aliquot of the effluent was sampled to measure HTO activities with a liquid scintillation counter (Packard 3170 AB/TR, Shelton, CT, USA). During the first upflooding process, both control valves were open to reservoir No. 2, and the flow rates through the column were controlled at approximately 5.0 ± 0.1, 3.0 ± 0.1, and 1.0 ± 0.1 mL/min. After 3 PV of the liquid phase passed through the column in the upflooding process, the HTO was flushed out from the column in a downflooding process where switch valves were open to reservoir No. 1 (GW). The effluent was collected using an auto-fraction collector every 60, 100, and 300 s. Table 2 shows that this up/downflooding process was set up and repeated two more times. A 10 mL aliquot of the sample was mixed with 10 mL of scintillation cocktail in a 20 mL polyethylene counting vial. All chemicals used in the experiments were of analytical purity, and de-ionized water (DIW) was used throughout the experiments.  The transport model was calibrated/validated by designing ADE experiments at different flow rates [33]. I − and IO 3 − , as well as disodium potassium iodide (KI) and potassium iodate (KIO 3 ), were injected into the column as stable isotope tracers after the HTO in the column was completely washed out. The corresponding concentrations (C(t)) of I − and IO 3 − were determined via ion chromatography (IS-600, Thermo Fisher Scientific Inc., Waltham, MA, USA) for I speciation (I − or IO 3 − ), and the iodine concentrations in the solution were measured with an induced coupled plasma optic emission spectrometer (iCAP 7000, Thermo Fisher Scientific Inc., Waltham, MA, USA). Table 2 lists the results of a multi-stage ADE column test, which included a series of upflooding/downflooding processes, for I − and IO 3 − . During the advection-dispersion experiments (ADEx) period for I − and IO 3 − , the frequency of collection depended on the flow rate, that is, every 60 s for approximately 3 mL/min. The corresponding concentrations (C(t)) of I − and IO 3 − in the effluent were then determined, and the respective breakthroughs were plotted for analysis. To ensure that the ADEx was finished, the upflooding process was stopped when the relative concentration (C/C 0 ) reached approximately 1.0, and the downflooding process was stopped when C/C 0 decreased to 0.05-0.1.

Mathematical Model and Parameter Estimations
STANMOD is a Windows-based computer software package for evaluating solute transport in porous media using analytical solutions of the convection-dispersion solute transport equation. The software package includes a modified and updated version of the CXTFIT code for estimating solute transport parameters using a nonlinear least-squares parameter optimization method. Residuals are a measure of the degree to which experimental results (Ce) deviate from the predicted values (Cp), N is the number of experimental data, and the root mean square error (RMSE) is defined as The HYDRUS-1D program numerically solves the Richards equation for variably saturated water flow and advection-dispersion type equations for heat and solute transport. It uses a numerical solution for numerical fitting. The transport equations include provisions for nonlinear nonequilibrium reactions between the solid and liquid phases. In addition, physical nonequilibrium solute transport can be accounted for by assuming a two-region, dual-porosity type formulation that partitions the liquid phase into mobile and immobile regions. As part of the inverse solution, HYDRUS produces a correlation matrix that specifies the degree of correlation between the fitted coefficients. An important measure of the goodness of fit is the r 2 value for regression of the observedŷ i versus fitted y i values: The r 2 value is a measure of the relative magnitude of the total sum of squares associated with the fitted equation; a value of 1 indicates a perfect correlation between the fitted and observed values. where w i is weighting factor for the two overlapping regions of the model,ŷ i is fitting value, y i is observed value.

Elemental Analysis by SEM-EDS
XRD, micro-X-ray computed tomography, polar microscopy, and analysis in previous works [31,32,34] showed that granite contains crystalline particles because it is igneous rock. In addition, we obtained and compared the corresponding images of SEM-EDS and elemental mapping analysis ( Figure 2). Comparison of the mapping area and mineral components in granite showed that O, Si, Al, Mg Fe, K, Na, Ti, Mn, and Ca are the major mineral components, which agrees with previous studies that used different analysis methods, such as Rutherford backscattering spectrometry [34].

Experimental BTCs of HTO, I − , and IO 3
− Before using I − and IO 3 − , a nonreactive radiotracer (HTO) was applied to characterize the major physical transport processes in the proposed dynamic column system with a PV (1 PV = 28 mL). Figure 3 shows that the HTO, I − , and IO 3 − BTCs at different flow rates reached 1 (C/C 0 = 1) from a series of sorption/desorption (up/downflooding) processes. Our experimental results indicated the absence of dead-end or dead pores to block the HTO, I − , and IO 3 − pathways, and water was saturated in each pore space in the compacted granite powder. It can be seen from the figure that since the injection method of solute is to continue to inject for a period of time and then stop the injection, all penetration curves are "s" type (divided into "s" type in the injection adsorption stage of the first half and "s" type in the desorption stage of the second half). Specifically, when the injection amount of solute reaches about one column pore volume, the ratio of solute activity concentration in the effluent quickly reaches the peak value, and then the concentration rapidly decreases to zero after continuous injection of clean water.

Fitting BTCs of Nonreactive Tracer (HTO)
The nonreactive tracer (HTO) was designed and applied in multi-stage ADEx to identify the reliability of our ADE device to build adequate reaction models through a calibration/validation process. Figure 4 (Fit-1 to -3) and Table 3 show that the equilibrium model in STANMOD-Fit (S) and HYDRUS-1D-Fit (H) were applied to fit the experimental data to obtain the parameters for HTO transport in the granite. The parameters (Fit-1 to -3) obtained by the two fitting methods showed similar results, indicating that the filling granite in the column showed good homogeneity. In addition, the highly reliable and accurate experimental data for HTO were identified according to the dispersivity (D L = 0.270-0.289 cm), retardation factor R (=1), and distribution coefficient K d (=0) in our ADE system by fitting with STANMOD and HYDRUS-1D.   Figure 5 shows the BTCs of anionic I − and IO 3 − in the granite column at different flow rates, and only the BTCs of IO 3 − (Fit-4 to -6) were similarly displayed and compared with the BTC of HTO. When the elution concentration ratio (C/C 0 ) of IO 3 − reached approximately 0.5 at 1 PV, R was almost equal to 1 and similar to HTO. Table 4 shows the fitting parameters of IO 3 − . The R (STANMOD) of IO 3 − and distribution coefficient K d (HYDRUS-1D) were almost equal 0 (=0), showing the non-reactive or adsorptive behavior of IO 3 − in granite. As shown in Figure 5d and Table 5, the BTC C/C 0 was less than 1 PV (Fit-7 and -8) of I − , and the R of I − was lower than 1. An obvious anion repulsive effect occurred, which accelerated the I − transport compared with HTO and IO 3 − transport. In fact, for the BTC of I − , when the equilibrium model was used for fitting, when Q = 5.0 ± 0.1 mL/min, the analytical solution could be used for fitting to obtain a good effect. When the numerical solution was used for fitting, the fitting effect was poor. This phenomenon did not occur when fitting the penetration curve with Q = 1.0 ± 0.1 and 3.0 ± 0.1 mL/min. In terms of I − fitting parameters, the R of I − obtained by fitting 5.0 ± 0.1 mL/min (Fit-7) and Q = 3.0 ± 0.1 (Fit-8) curves with analytical solution was less than 1, possibly indicating obvious anion exclusion. After fitting with the numerical solution in Table 5, the K d of the curve with Q = 5.0 ± 0.1 mL/min was greater than 1 by default (Fit-7). Thus, the fitting correlation was not good enough (r 2 = 0.77) compared with Fit-8 and -9 (r 2 > 0.95). The K d of I − Q = 3.0 ± 0.1 mL/min was smaller than that of HTO, and I − transport in granite could not be explained by fitting the equilibrium models (STANMOD and HYDRUS−1D). Table 6, Fit-10 and 11 (r 2 > 0.95) show that the two-region nonequilibrium model was developed [33] and applied in this work to fit the BTC curves of Q = 5.0 ± 0.1 and 3.0 ± 0.1 mL/min (STANMOD) and Q = 5.0 ± 0.1 mL/min (HYDRUS-1D) and understand the I − transport in granite. The fitting parameters showed that the water contents of immobile liquid phases (θ im ) increased with increasing flow velocity (Fit-10), indicating that the anion exclusion effect was more obvious. In addition, the θ im was larger than the analytical solution possibly because of the numerical algorithm.
The migration behavior of HTO and I − in the same medium was compared microscopically to explain the influence of anion exclusion of I − . Figure 6a describes the migration of HTO in the groundwater environment. As shown in the figure, HTO can pass through all water channels under the action of convection dispersion and finally fill the whole medium. Figure 6b shows the migration of I − in the groundwater environment. According to the "electric DDL theory", the anion concentration in a single pore increases exponentially with the distance from the pore wall. When I − migrates in a wide water passage, it will converge to the center due to anion exclusion and will not pass through a narrow passage, forming "immobile regions". Finally, I − will select the channel with a wide water surface to migrate, resulting in acceleration.   Figure 7 shows a comparison of the migration behaviors of I − and IO 3 − under different flow rates. According to the theory of the "electric double layer model", a negative ion layer forms when the dielectric surface is negatively charged, and a positive ion layer forms outside due to electrical attraction. The closer it is to the dielectric surface, the more negative its potential becomes (the more obvious the repulsion). Comparison of the ion radii of I − and IO 3 − showed that the ion radius of I − (3.31 Å) was smaller than that of IO 3 − (3.74 Å) [35]. Therefore, I − can pass through the diffusion layer and become closer to the medium surface during migration. At the same time, the greater the migration speed of ions, the greater the kinetic energy generated, the closer the energy to the medium surface, the stronger the anion exclusion effect, and the more obvious the acceleration during migration. However, it is difficult for IO 3 − to become close to the medium surface because of its large ion radius and weak anion exclusion. Thus, the acceleration during migration was not obvious.

Conclusions
The transport behavior of HTO, I − , and IO 3 − in granite was studied by ADE column experiments. Comparison of the BTCs of HTO, I − , and IO 3 − under different flow rates showed that the anion exclusion effect occurred in I − at Q = 5.0 ± 0.1 mL/min during fitting STANMOD and HYDRUS-1D when R < 1. The key parameters of HTO, I − , and IO 3 − were obtained by fitting the equilibrium and nonequilibrium models of the analytical method (STANMOD) and numerical inversion (HYDRUS-1D). The conclusions are as follows: 1.
The BTCs for HTO were symmetrical at various flow rates (1.0 ± 0.1-5.0 ± 0.1 mL/min), and no significant differences were found in the accessible porosity and dispersivity values of HTO. Therefore, the experimental apparatus was highly reliable, and the granite samples filled and compacted in the column were almost homogeneous.

2.
Comparison with the BTCs of HTO, I − , and IO 3 − at different flow rates showed that obvious anion exclusion only occurred in I − transport by increasing the flow rate from 1 to 5. In fact, the anion exclusion (R < 1) only occurred for I − at a flow rate of 5 mL/min, and a relative Coulomb's repulsive force may be caused by the smaller hydration radius of I − (3.31 Å) than that of IO 3 − (3.74 Å) according to electric DDL theory. 3.
The equilibrium and nonequilibrium transport models were used and compared to identify the mobile/immobile zones in the compacted granite column. The anion exclusion effect depended on the immobile zones in the column. In sum, the nonequilibrium model can well characterize the immobile regions of anion exclusion, which is obviously proportional to the value of immobile regions. 4.
Since the experiment was conducted in laboratories, the environmental conditions were quite different from the actual environment. Whether the research results can be applied to the actual plant site needs to be demonstrated. It is recommended to carry out field experiments as soon as possible.