Analysis of Uranium Sorption in a Laboratory Column Experiment Using a Reactive Transport and Surface Complexation Model

The transport and retardation of radioactive elements in hyper-alkaline conditions of radioactive waste repositories is a challenging field that is still poorly understood. In this study, the transport and attenuation of uranium in a column experiment was modelled by considering kinetic reactions, advection–dispersion and chemical/physical retardation processes. The modelling was first performed for three alluvium samples from Yucca Mountain in circumneutral pH to moderately alkaline conditions. Sorption of uranyl (UO22+(UVI)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{UO}}_{2}^{2+} ({U}_{\mathrm{VI}})$$\end{document}) was found to strongly depend on the surface complexation model assumed, with no significant removal of UVI\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${U}_{\mathrm{VI}}$$\end{document} by precipitation or ion exchange process. The surface/edge site reaction of Al-hydroxyl group in kaolinite was shown to have a high affinity for uranyl adsorption, while the hydrous ferric oxide edge on hematite adsorbed most of the uranyl ions. The model was then used to interpret uranium transport in a laboratory column filled with Hollington sandstone under hyper-alkaline (pH 13) conditions. The simulation results show that uranium adsorption on the Al-hydroxyl edge of kaolinite exceeds adsorption by the calcium silicate hydrate phase. This result may reflect the lack of surface complexation parameters for calcium silicate hydrate minerals. Hence, further studies are required in the field of surface complexation reactions for calcium silicate hydrate phases. Sorption of uranyl (UO22+(UVI)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{UO}}_{2}^{2+} ({U}_{\mathrm{VI}})$$\end{document}) was found to strongly depend on the surface complexation model, with no significant removal of U_VI by precipitation or ion exchange process. The aluminol surface edges in kaolinite were shown to have a higher affinity for uranyl adsorption, while the hydrous ferric oxide edge on hematite adsorbed most of the uranyl ions. Uranium adsorption on the aluminol edge of kaolinite exceeds adsorption by the C-S-H phase. This result may reflect the lack of surface complexation parameters for C-S-H minerals. Sorption of uranyl (UO22+(UVI)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{UO}}_{2}^{2+} ({U}_{\mathrm{VI}})$$\end{document}) was found to strongly depend on the surface complexation model, with no significant removal of U_VI by precipitation or ion exchange process. The aluminol surface edges in kaolinite were shown to have a higher affinity for uranyl adsorption, while the hydrous ferric oxide edge on hematite adsorbed most of the uranyl ions. Uranium adsorption on the aluminol edge of kaolinite exceeds adsorption by the C-S-H phase. This result may reflect the lack of surface complexation parameters for C-S-H minerals.


Introduction
Uranium is a radioactive and toxic element which, besides its natural occurrence, can be distributed to the environment through human activities such as oil and gas production, mining processes, and the nuclear industry (Campos et al. 2011;Chandrajith et al. 2010;Hasche-Berger 2006, 2008;Ricka et al. 2010). Those industrial processes can enrich the radionuclides activities of uranium and its decay products known as NORM (naturally occurring radioactive material) (Protection 2003). In an oxidising environment the most stable and soluble form of uranium is (Uranyl UO 2+ 2 (U VI ) ), which can react with sulphates, carbonates and nitrates to form complex phases with altered speciation and transport characteristics (Zielinski et al. 1997;Meinrath 1998;Tutu et al. 2009;Ribera et al. 1996;Grenthe et al. 1992). In subsurface applications the aqueous concentration and mobility of U VI ions without a complexing ligand are widely recognised to be controlled by sorption at mineral surfaces (Barnett et al. 2002;Prikryl et al. 2001;Turner et al. 1996a; Thomson et al. 1986;Payne and Waite 1991). The interaction between dissolved radionuclides and minerals in the host rock (sorption at solid/solution interface) is critical when evaluating the immobilisation of radionuclides.
Uranyl sorption on mineral surfaces can involve multiple processes (binding sites), such as absorption, adsorption, ion exchange and edge/surface complexation. These processes can significantly influence the transport of pollutants and radionuclides in soils and rocks. For example, adsorption refers to the adherence of ions to the solid surface, while absorption implies uptake of ions into the solid. Conversely, ion exchange involves the substitution of one ion for another over the solid surface (Appelo et al. 2005). In earlier studies, sorption models use the K d (distribution coefficient) approach, which sums all interactions between the solid/water interaction (Ticknor 1994;Missana et al. 2008;Glynn 2003;Davis et al. 2004;Curtis et al. 2004).
In recent years, a mechanistic approach has been developed based on several absorption/adsorption laboratory experiments conducted on site-specific clay materials to generate a surface complexation model (SCM) that can describe the migration behaviour, stability constants and stoichiometry of different metal/metalloid ions reactions. The SCM can describe and simulate sorption between aqueous uranyl and the surface of clay minerals (sorbent) by fitting the thermodynamic data from other referenced studies under different chemical conditions for various aqueous (uranium) species (Davis 2001;Davis et al. 2001;Bradbury and Baeyens 2005;Zachara and McKinley 1993;Kim 2001;Davis et al. 2002;Koretsky 2000). Through the mass action equation, the surface complexes formed can be then integrated within the reactive transport models (Curtis et al. 2006;Papini et al. 1999).
In high-heat-generating waste-disposal concepts, swelling clays (e.g. bentonite) will be widely used, in part, because of their capability to sorb heavy metals. These clays can sorb heavy metals due to their high osmotic swelling capacity and large specific surface area (Barnett et al. 2000(Barnett et al. , 2002McLing 1998). Moreover, clay minerals can simultaneously maintain a fixed, usually negative charge within the particle structure which encourages ion exchange within the interlayer's spaces, while the variable charge is generated at the clay edges, which may be positive or negative depending on the pH value. Therefore, both cations and/or anions can potentially be sorbed to neutralise the structural charge. Three common surface complexation models are usually used to fit the experimental results: the double/diffuse layer model (DLM), the triple-layer model (TLM), and the constant capacitance model (CCM) (Waite et al. 1994;Turner et al. 1998Turner et al. , 1996aHiemstra et al. 1989;Fletcher and Sposito 1989). Due to the low hydraulic conductivity of clay soils, most investigations have been carried out in static systems using batch experiments with powdered material left in contact for periods of time with cement leachate (Claret et al. 2002;Ramıŕez et al. 2002;Fernández et al. 2009;Takahashi et al. 2007). Fewer studies have performed dynamic experiments using diffusion set-ups or advective experiments (Adler et al. 1999).
The aim of this work is to establish a reactive transport model with a series of kinetic and equilibrium reactions and incorporate uranyl aqueous speciation and surface charge measurements that control uranium removal from the solution at room temperature and fixed pH value. The modelling results were obtained using the PHREEQC geochemical code. In past numerical modelling attempts, when analysing nuclear waste repositories, most authors assumed the host rock's key properties to be constant. Though, in this research, an adequate approach was used to study the effects of variable porosity, sorption, reactive surface area and pore volume on improving the modelling of rock alteration in the system.

Experiment
In this study, a reactive transport model combined with a diffuse layer model (DLM) was used to analyse the adsorption of uranyl ions from aqueous solution in column experiments. The objective was to determine the key parameters that control the sorption of this radionuclide at the solid-solution interface in the host rock of a candidate geological repository. The model was first applied to three alluvium samples from Yucca Mountain that are rich in smectite clay and equilibrated with circumneutral to moderately alkaline solutions. Once the model was validated against the experimental results, it was then used to interpret the transport of uranium in a highly alkaline solution through Hollington Sandstone, a material mainly composed of silicate minerals with minor amounts of kaolinite and hematite. Hollington Sandstone has been used for a mineral analogue for host rocks of UK ILW/ LLW disposal.

Yucca Mountain Alluvium
The Yucca Mountain project in Nevada has been proposed as a geological repository site for the disposal of radioactive wastes. A series of sorption studies (batch, column, and in situ field transport) have been conducted on the alluvium soil to demonstrate its capability to retard the transport of radionuclides in the subsurface.
Three flow-through column experiments were conducted in the early 2000's with continuous flow and fully saturated conditions using alluvium from the southern area of Yucca Mountain (Scism 2005). The columns were filled with alluvium that had been wet sieved to retain the 2000-75 µm size fraction. The experimental parameters for the three column tests are shown in Table 1. The key mineral phases in the alluvium were identified and quantified using X-ray diffraction (XRD) following the procedure published by Chipera and Bish (2002). Table 2 shows the mineralogy of the alluvium in each of the three columns. The background solutions (groundwater) used in the experiment were taken from water wells drilled close to each alluvium sample wells. The water chemistry and ion concentration data for each column are shown in Table 3. (A concentrated UO 2 (NO 3 ) 2 solution was diluted in the groundwaters to produce a uranium tracer solution with a concentration of 1 × 10 -6 M. The columns were initially allowed to equilibrate with groundwaters and then the uranium tracer solution was injected at an initial flow rate of 10 ml/h, which decreased to 3 ml/h as the experiment continued. The fluid sample at the column outlet was placed in a reciprocating shaker, followed by centrifugation to separate the solids. Simultaneously, control samples that contained only a tracer without alluvium were also shaken and centrifuged to estimate the amount of uranium sorbed to the tube walls. Aqueous uranium concentrations were determined by Liquid Scintillation Counting (LSC) (the concentration in the tracer solution was measured before and after the experiment). Further details of the experimental procedure and protocol are provided in (Scism 2005).

Hollington Sandstone
In this experiment, uranium in a high-alkaline pH solution will migrate through the sandstone sample and alter the solution composition due to the dissolution of primary silicate minerals in the rock. Later, the solution will be equilibrated by the precipitation of calcium silicate hydrate (C-S-H) phases or C-(Al)-(K)-S-H phases whenever Al and K ions are present in the solution.
Hollington sandstone is part of the Lower Triassic Bromsgrove Sandstone Group. It is mainly composed of silicate minerals with minor quantities of clay and hematite. Its porous nature allows the study of rock-cement-leachate reactive transport reactions in strata with a similar mineralogical composition to that a proposed geological radioactive waste repository. Column experiments on crushed Hollington sandstone permeated with a synthetic young cement leachate at a temperature of 50 °C (Small et al. 2016) are used to predict the transport of aqueous U VI in Hollington sandstone. In this experiment, the aim was to investigate the progression of uranium sorption on the surfaces of the C-S-H phases in an alkaline environment. In fact Langmuir (1997) concluded that hyper-alkalinity conditions can encourage the sorption of uranyl ( UO 2+ 2 ) cations due to the availability of neutral and negatively charged surfaces. Furthermore, the C-S-H phases have a high surface area and high retention capacity for radionuclide migration, especially for uranyl (VI) cations due to poor crystallisation (Johnson 2004;Gougar et al. 1996;Ma et al. 2019). Meanwhile Tits et al. (2011) found that C-S-H phases can absorb uranyl cations by forming inner-sphere surfaces (surface complexation) with the silonal edges. Tables 4 and 5 show the mineralogical composition of the sandstone sample used in the experiment and the chemical composition of the synthetic young cement leachate, respectively. Further details of the experimental set-up and procedures are provided in Baqer et al. (2021).

Modelling Approach
The geochemical speciation code PHREEQC (version 3.6.1) (Parkhurst and Appelo 2013) was used to simulate the transport and sorption of uranium in the column experiments. PHREEQC is a geochemical software developed by the United States Geological Survey. It includes multiple databases users can use or modify and apply to various aqueous models, including geochemical and industrial applications. Adding to its geochemical modelling tool, it can simulate 1D transport in porous media, making it applicable for the resolution of this research. The software comprises a vast database and geochemical equations describing the equilibria or kinetic interactions of aqueous solutions with solids, gases, and minerals, with the advantage of easily adding or modifying necessary reactions if needed. In PHREEQC, chemical reactions are integrated based on a built-in or user-defined rate expression, with the possibility of simulating time-dependent geochemical reactions that depend on solute temperature, pH and saturation ratios of minerals. Regarding kinetics reactions, the built-in interpreter contains rate expression (which can also be user-defined) with a Runge-Kutta scheme that solves the differential equations and change the solution speciation within a specific time interval. In addition, the user can define an error tolerance which allows the software to reduce the time interval of the simulation. This scheme effectively achieves equilibrium when multiple kinetic reaction rates are defined, and the rate of each reactant is changed during the reaction.
In this work, the kinetics of dissolution and precipitation, equilibrium reaction and porosity evolution have been incorporated in the coding. A hybrid mixed kinetic-equilibrium approach was used to overcome the shortage and uncertainties of some kinetic reaction parameters (e.g. reactive surface area, specific dissolution/precipitation kinetics) (Chen and Thornton 2018;Van der Lee 1998;Bethke 1996). The LLNL thermochemical database was used as a starting point to compile the details of different chemical reactions (aqueous reactions, mineral dissolution/precipitation, surface complexation, ion exchange) (Delany and Lundeen 1990). This database seemed to be the best option available since it has the thermodynamic information (reactions and equilibrium constants) for variety of minerals and aqueous species, especially carbonate minerals that have major influence on key kinetic reactions and putative reactive pathways that controlled the primary mineral dissolution, mineral evolution, secondary mineral formation and column effluent evolution. Uranium thermodynamic data were obtained from the "Second update on the chemical thermodynamics of uranium, neptunium, plutonium, americium and technetium" published by the OECD Nuclear Energy Agency (Grenthe et al. 2020). Surface complexation and ion exchange reaction were obtained from the literature and then added manually to the database.

Equilibrium/Kinetic Modelling of Dissolution and Precipitation
The sorption of a soluble radionuclide can be highly affected by dissolution/precipitation processes which occur at mineral surfaces. In some circumstances, the dissolution of some mineral phases can result in the nucleation of clay minerals or hydrous oxides as a coating on the mineral surfaces, which will increase the number and density of sorption sites (Ticknor 1994). Also, precipitation process can sometimes involve radionuclides, removing them from the solution by formation of a secondary phase. The alluvium soil sample contains 11 minerals with different percentages in each column. To ensure consistency in the results and to highlight the sorption effect of each mineral, all the mineral phases were included in the modelling except for "Opal-CT" since it is only found in column 1. The higher percentage minerals (quartz, plagioclase, k-feldspar, clinoptilolite, cristobalite, and smectite) were modelled kinetically while the rest were included in the equilibrium reaction. Equation 1 was used to calculate the rate of mineral dissolution/precipitation reactions (Appelo et al. 2005): where R is the overall dissolution/precipitation rate (mol L −1 s −1 ), k is the specific dissolution/precipitation rate (mol m −2 s −1 ), A is the effective surface area of the mineral (m 2 g −1 ), V is the pore volume (L), M is the moles of solid at a given time, M 0 is the initial moles of solid, and (IAP∕K) is equal to the saturation ratio (SR) value of the mineral where IAP is the Ion Activity Product and K is the equilibrium constant. In the above equation, the  (Carroll and Walther 1990) molar concentration (M∕M 0 ) has n power, which is a function of the initial crystal grain size distribution that affects the solid dissolution rate. A value of 2/3 is assumed in this work (this is a typical value for uniformly dissolving cubes or spheres in the monodisperse population) (Appelo et al. 2005). Table 6 shows the specific kinetic rates and the reactive surface areas for the kinetically modelled minerals. To investigate the effect of precipitation on uranium removal from solution, four uranium mineral phases were added in the equilibrium reactions (In PHREEQC) during the simulation of the Yucca Mountain and Hollington sandstone experiments (CaUO 4 , Uranophane, Becquerelite, Rutherfordine). These phases were chosen based on the common stable uranium minerals found in the literature (Atkins et al. 1988;Wronkiewicz et al. 1992;Gorman-Lewis et al. 2008;Felipe-Sotelo et al. 2017). Also, Tobermorite-14A, CSH-gel, CaUO 4 , and Na 2 U 2 O 7 were included in the modelling of the Hollington sandstone experiment as a possible precipitating phase in a high-alkaline environment.

Dynamic Porosity and Reactive Surface Area
In the model, a dynamic porosity was used to emulate changes in the porosity that result from dissolution/precipitation reactions. Equation 2 shows the correlation between mineral mass and surface area during flow of the tracer solution in the column.
where m t is the total mass of the rock sample, A 0 is the initial surface area ( m 2 g −1 ), and X i is the mass fraction of each mineral in the rock sample (Beckingham et al. 2016). The surface area is then coupled with the evolved porosity through Eq. (3) (Lichtner 1988): where A t=0 r is the reactive surface area of the mineral at the initial porosity ( t=0 ) , and the mineral kinetic rate can be calculated using Eq. (4). The changed porosity can then again be calculated from the volume of the dissolved/precipitated moles of mineral. Further details concerning the mathematical derivation of the dynamic porosity model are provided in Baqer et al. (2021).

Transport Process
The transport of the uranium tracer solution through the column was modelled using the one-dimensional (1D) mixed cells concept with flux type boundary conditions and a mass entering the column per unit time. This is applicable to a laboratory column with a diameter much smaller than the column length (Table 1) because of the near-zero concentration gradient at the column end (Appelo et al. 2005). The column was divided into ten equal length cells along the flow path, with additional cells at each end for inflow and outflow. The simulation was performed for a series of time steps in which each step is equal to the time taken for the pore solution in one cell to move into the next cell. The simulation was performed for 700 h, and the advection-dispersion process through the column was modelled using the equation of advection-reaction-dispersion described below (Eq. 5) (Nardi et al. 2014): where C i is the total concentration of component i in water (mol Kg w −1 ), t is time (s), x is the distance (m), D L is hydrodynamic dispersion coefficient (m 2 s −1 ), im is the stoichiometric coefficient of component i in kinetic reactant m (dimensionless), R m is the overall dissolution/precipitation rate of kinetic reactant m (mol Kg w −1 s −1 ), and n is the number of kinetic reactants. For each time step, advection was modelled first, followed by dispersion, then finally the chemical reactions were modelled, so the solution composition in each cell could be updated. The reader is referred to Baqer et al. (2021) for more details concerning the transport model. In the Yucca Mountain experiment, different flow rates for the pulse solution were used for each column. Based on the time step in each column, the velocity will vary, and therefore each column will have a different dispersivity value. Table 7 shows the modelling input parameters for the transport process simulation in the three columns.

Uranium Speciation
The aqueous speciation of the uranyl ( UO 2+ 2 ) ion is strongly influenced by pH, as it determines the stability of different ion complexes and the distribution of the surface sites, which will consequently shape the sorption mechanism . For example, where CO 2 is present, neutral uranyl-carbonate species start to form at pH values > 6.5, and negatively charged uranyl-carbonate species dominate in alkaline conditions (Nair et al. 2014). Formation of negatively charged complexes will affect adsorption to partially ionised surface/edges sites (Guimarães et al. 2016;Ticknor 1994). Table 8 lists equilibrium constants for uranium aqueous speciation reactions at 25 °C and pCO 2 = 10 -3.5 atm (Grenthe et al. 2020).

Sorption
Smectite contains several types of sorption site with high affinity for aqueous uranium complexes (Barnett et al. 2000;Davis et al. 2004), and kaolinite contains aluminol (5) 0.3 × 10 -9 0.3 × 10 -9 0.3 × 10 -9 surface/edges that also have high affinity for uranium species (Borovec 1981;Kohler et al. 1992). Furthermore, the presence of iron ions in hematite can further enhance uranium sorption, hence it has also been included in the surface complexation model Dong and Wan 2014). Most published literature model sorption with cation exchange and surface complexation mechanisms (Turner et al. 1996a;Missana et al. 2008;Bachmaf and Merkel 2011;Korichi and Bensmaili 2009;Nair et al. 2014). The advance in this study is that it models both ion exchange and surface complexation as multi-site complexation phenomena at the clay-water interface (Fig. 1).

Ion Exchange
Uranyl sorption at smectite interlayer cation exchange sites is included in the model. It assumed that these sites are initially occupied by either Ca 2+ or Na + cations [Na + 26.1 Fig. 1 Multi-site complexation model of uranium on smectite clay is readily displaced by UO 2+ 2 , whereas Ca 2+ has a similar affinity to interlayer cation exchange sites to uranyl (Tsunashima et al. 1981;McKinley et al. 1995;Hiemstra et al. 1989)]. The stoichiometry and equilibrium constants assumed for these cation exchange reactions are shown in Table 9. The site density (X − ) is assumed to equal the cation exchange capacity (CEC), which for smectite equals 810 meq/kg (Appelo et al. 2005).  2014) is used for uranium sorption onto various surface and edge sites on clay minerals to overcome the limitation of reactive solute transport using the average K d value (laboratory calculated) approach (which depends on, and is limited to, a specific water composition) (Glynn 2003;Reardon 1981;Zhu 2003;Bethke and Brady 2000;Davis et al. 1978;Kent et al. 1988). The model assumes that the mineral surface includes a group of functional hydroxyl surfaces denoted as (≡ SOH) , and that the sorption of uranium is highly dependent on the behaviour of these functional groups at surface and edge sites on clay minerals. Such sites are amphoteric, so protonation and deprotonation equilibrium reactions of these sites are included in the model (Missana et al. 2008) where K 1 and K 2 are the equilibrium constants (Table 10).

Surface Complexation
In this study a double/diffuse layer SCM model was used, and the moles of each aqueous species are computed based on a constant thickness for the diffuse layer. Before applying the DLM for the adsorption of uranyl ions onto clay surfaces, several variables must be defined: (1) the chemical reaction at the clay-water interface, (2) the surface complexation equilibrium constant (Log k), and (3) surface site density and amount of available binding sites.
Usually, for smectite, the binding of uranyl aqueous complexes takes place on the aluminol and silanol edges, and on the aluminol edge site for kaolinite (Bachmaf and Merkel 2011;McKinley et al. 1995;Turner et al. 1996b;Zachara and McKinley 1993). The amount of edge/surfaces associated with surface complexation is usually 10-20% of the sorption sites (Anderson and Sposito 1991;Pabalan and Turner 1996;McKinley et al. 1995;Hiemstra et al. 1989). Table 10 shows the surface edge site reactions on aluminol/silanol hydroxyl groups included in the PHREEQC modelling to best fit the experimental uranyl sorption results. Dzombach and Morel (1990) recommend using 10% of the total specific surface area measured by the N 2 -BET method to account for (6) ≡ SOH + 2 ⇔≡ SOH + H + K 1 (7) ≡ SOH ⇔≡ SO − + H + K 2 UO 2+ 2 + CaX 2 = UO 2 X 2 + Ca 2+ 0.049 Guimarães et al. (2016) the area of the crystallite surface edge sites, or site density value of 2.3 sites/nm 2 for all minerals. Therefore, the total number of moles of binding site surfaces was calculated from the weight percentage and surface area of each mineral in the bulk sample of alluvium in the column. Note that the surface complexation reaction of carbonate complexes has been ignored in this study since they have high solubility.

Yucca Mountain Alluvium
The final composition of the effluent solution was calculated after a series of geochemical reactions combined with fluid flow and uranium sorption. The results were obtained by implementing multiple binding site ion complexation containing fixed and edge aluminol/ silonal sites. Table 11 shows the predicted composition of the site groundwater after equilibration with mineral phases in Yucca Mountain alluvium and atmospheric carbon dioxide. In all cases the initial breakthrough of uranium almost occurs at the same time as the HTO-tracer (Fig. 2). The modelling results describe the adsorption of uranyl well: the uranium concentration starts to increase in the existing fluid until it reaches a maximum concentration (C/C 0 ) of 0.0078 mol L −1 , 0.023 mol L −1 , and 0.08 mol L −1 for columns 1, 2 and 3, respectively (Fig. 3). However, the highest value of the uranium concentration was still much lower than the pulse solution concentration (0.24 mol L −1 ), implying that uranium retardation happens simultaneously with the fluid flow. The simulated saturation indices of uranium phases (CaUO 4 , uranophane, becquerelite, rutherfordine, Fig. 4) indicate that they are undersaturated and the precipitation of solid uranium phases is unlikely to occur in the experiment. Therefore, precipitation can be excluded as a possible uranium removal mechanism in this experiment. This is usually the case since the nucleation kinetics of solid phases is considered to be very slow to affect uranium removal significantly (Dangelmayr et al. 2017). Consequently, the decrease in uranium concentration is due to the sorption of uranyl ions by ion exchange and surface complexation reactions only. Figure 5 shows the number of moles of uranium adsorbed by cation exchange processes. As the experiment proceeds, the pulse solution is exposed more to the smectite mineral; thus, more uranyl is exchanged on the smectite binding sites. However, the amount of adsorbed uranyl is minimal throughout the experiment, which agrees well with other studies that highlight the dominance of ion exchange in the region of low pH only (McKinley et al. 1995;Zachara and McKinley 1993;Davis and Kent 2018). Therefore, uranium binding to a fixed site by ion exchange process is most likely not the dominant sorption mechanism in this model as well. Figure 6 illustrates the results of the four surface edges sites ( AlOH Smectite , SiOH , AlOH Koalinite , FeOH ), which confirm the modelling approach for uranium sorption of binding to a variable charge site by surface complexation. As expected, uranyl ions form an inner surface complex on smectite, kaolinite and hydrous ferric oxide. When assessing the adsorption results of each surface site, column 3 shows the highest adsorption value in all four sites, followed by column 2 and finally column 1. This agrees well with the highest smectite and hematite weight per cent in column 3 and the lowest weight percentage in column 1. Moreover, from Table 11, column 3 has the highest ionic strength while column 1 has the lowest, which agrees with the results of other literature that shows a more elevated surface complexation potential with higher ionic strength (Chisholm-Brause et al. 2001;Hennig et al. 2002;Sylwester et al. 2000;Bauer et al. 2001;Korichi and Bensmaili 2009;Schindler et al. 2015;Turner and Sassman 1996).
The number of adsorbed moles on the aluminol edge (≡ AlO − UO + 2 ) is much higher (6 orders of magnitude) than the number of adsorbed moles on the silonal sites ( ≡ SiO − UO + 2 ) , which agrees well with higher uranyl affinity towards aluminol edges. This may be justified by the lower tendency of the ( ≡ SiO − UO + 2 ) group to donate its oxygen and form inner-sphere surface species Kowal-Fouchard et al. 2004). Even with the high weight per cent of smectite clay in the alluvium samples, it has the lowest sorption capacity for uranium. The aluminol surface edge sites in kaolinite have shown a higher affinity for uranyl adsorption, which is also reported in other studies (Borovec 1981;Kohler et al. 1992;Payne et al. 2004) for a pH range between 5-9. This can be related to the higher number of the exposed surface sites at the Al-octahedral sheet (Kaolinite has a 1:1 clay structure), which results in greater uptake of uranyl ions. In all three columns and among the four surface complexes, the ≡ FeOH site has the highest number of adsorbed uranyl moles. The same behaviour has also been reported in other studies (Dzombak and Morel 1991;Waite et al. 1994;Sherman et al. 2008;Hsi and Langmuir 1985), in which for the region of neutral and alkaline pH, a very small mass weight percentage (≈ 1%) of iron hydroxides (goethite, hematite, ferrihydrite) is a major sink for uranium. More importantly, (Liger et al. 1999) found that in the presence of hematite, uranium reduction through sorption can occur within hours at neutral pH; the reaction kinetics between the adsorbed uranium and ferrous iron is enhanced by hematite according to a first-order pseudo-kinetic law.
The distribution of uranium species in PHREEQC also shows the presence and dominance of Ca 2 UO 2 (CO 3 ) 3 , CaUO 2 (CO 3 ) −2 3 complexes. This is consistent with the view of several authors that neutral and negatively charged uranium-carbonate ions are unlikely to bond with the negatively ferric-hydroxide surfaces (Morrison et al. 1995;Geipel et al. 1998;Ho and Miller 1986). Moreover Fox et al. (2006), Dong and Brooks (2008) ,   Fig. 3 Experimental and modelling results for uranium breakthrough curves in Yucca Mountain alluvium columns 1, 2 and 3

Fig. 4 Saturation indices for uranium secondary phases in Yucca
Mountain alluvium columns 1, 2 and 3 1 3 Stewart et al. (2010), Nair and Merkel (2011a), Nair and Merkel (2011b) found that the availability of Ca and Mg ions in the solution can also shift the aqueous speciation of uranium towards more stable (neutral) or negatively charged Ca-Mg-ternary complexes. The simulation results also indicate a general dissolution of smectite, clinoptilolite, and cristobalite in all three columns (Fig. 7), which will result in an increase in the pore volume of each column and hence the porosity (Fig. 8).

Hollington Sandstone
The experiment was conducted in a CO 2 free environment which is the usual case in a cementitious geological repository due to the interaction with cement-Ca and precipitation of calcite (CaCO 3 ) (Disposal 2010;Vines and Lever 2013;Auroy et al. 2013). Therefore, in an oxidising environment, U VI ions would be expected to be in the form of uranyl hydroxide complexes with limited solubility compared with the carbonate species (Tits et al. 2011;Bourdon et al. 2003). The model simulations show that there is no significant precipitation of C-S-H or uranium composed phases (Fig. 9). Tobermorite-14A, CSH-gel, and CaUO 4 are all in multiple cycles of precipitation and dissolution that cancel out the effect of both processes. Meanwhile, Na 2 U 2 O 7 approaches but never reaches equilibrium. Therefore, no significant retardation is expected by the edge surfaces of the C-S-H phases or precipitation of uranium. This result is also represented in Fig. 10, which shows that in the case of no sorption, the uranium concentration at the column end reaches its initial injected concentration (0.24 mol L −1 ) and achieves a steady-state value. Conversely, in the full sorption model uranium breakthrough occurs simultaneously as in the previous case but with significant retardation of uranium.
As there was no significant precipitation of secondary C-S-H phases, an assumption was made to treat the smectite clay in the Hollington sandstone as a secondary C-S-H mineral. This assumption has also been made due to the lack of surface complexation parameters for C-S-H minerals. The modelled results for this assumption are shown in Fig. 11. Again, the results show that aluminol surface complexation on kaolinite dominates uranyl adsorption. So, even in high-alkaline conditions, it still overcomes the adsorption of both aluminol and silonal edges on the C-S-H phase. But it is worth mentioning that (Harfouche et al. 2006;Evans 2008;Pointeau et al. 2004) have experimentally found that uranyl (UO 2+ 2 ) can substitute for Ca 2+ in the interlayer of C-S-H phases. Unfortunately, this reaction could not be modelled in this study due to a lack of the relevant surface complexation and ion exchange parameters.  In the study of Korichi and Bensmaili (2009), the authors also found that increasing the uranium concentration in the pulse solution will decrease the adsorption percentage, and in the case of low initial concentration, uranyl ions will have high mobility at the beginning before adsorption occurs. This behaviour was perfectly captured in the modelling results of this study, where a small peak in uranium concentration occurs at the beginning of the experiment before it drops down once the adsorption reaction dominates (Figs. 10, 11). The applied uranium transport and sorption methods in this study provide insight into the key parameters that have the most significant impact on minimising uranyl mobility in a geological repository. It also indicates the required data needed to construct a consistent and reliable model for the reactive transport and sorption process of radionuclide migration in the geochemical application. Finally, the result of this study shows that the prediction of uranium migration is highly site-specific as it depends on the mineralogy and the geochemistry of the geosphere.

Conclusion
Modelling uranium retardation is still a challenging task with all the variables included in the process, such as characterisation of the sorbent material (mineral type and composition, density of sorption sites, surface area), solid-solution ratio, solution composition, carbon dioxide presence, and pH value. Usually, in the region of low pH and ionic strength, uranium is retarded by forming an outer sphere complex with the fixed charge surface (ion exchange process). In contrast, at neutral pH and high ionic strength, the formation of inner-sphere complexation due to amorphic surface/edge sites occurs (surface complexation) controls the sorption process. Precipitation of uranium phases is most likely to be neglected in modelling experimental uranium sorption since the nucleation kinetics of solid phases is considered very slow to affect uranium reduction significantly. In this study, reactive transport with a double-layer surface complexation model has been implemented to model uranium sorption in a column experiment with a dynamic transport process. The model was first applied to an alluvium soil sample from Yucca Mountain in a neutral pH environment. The result shows that the aluminol surface edge sites in kaolinite have a higher affinity for uranyl adsorption than both aluminol and silonal edges in smectite clay. Meanwhile, hydrous ferric oxide edge (≡ FeOH) on hematite adsorbed most of the uranyl ions. Subsequently, the model was applied to Hollington sandstone in a high-pH environment. The modelling simulation shows that there is no significant precipitation of C-S-H or uranium composed phases. Tobermorite-14A, CSH-gel, and CaUO 4 are all in multiple cycles of precipitation and dissolution that cancel out the effect of both processes. Meanwhile, Na 2 U 2 O 7 is almost close to equilibrium but never reached. Like Yucca Mountain, the aluminol surface edge sites in kaolinite overcome the adsorption of both aluminol and silonal edges on the C-S-H phase. It is worth mentioning that the surface complexation reaction on C-S-H phases was not modelled adequately due to the lack of modelling parameters and further studies are needed in this field.

Limitations of the mathematical model
Whether an experiment can be modelled correctly by numerical methods depends on (1) the experimental data obtained, (2) the modelling tools and methods, (3) the theoretical analysis, and (4) the relationship between the experiment and modelling (some parameters might be important for modelling, however, is not necessary for the experiment). One of the challenges of modelling this kind of experiment was the unknown parameters caused by the different scientific focus for the experiment. For example, one of the major uncertainties in modelling evolving porous media arises from the uncertainty in identifying the value of the reactive surface area, which has a first-order effect on the timing and the level of dissolution/precipitation reactions. Commonly, the rate of dissolution in the experiment is 1-3 orders of magnitude higher than in the natural system. This can be explained by the difference in reactive and total surface area between natural and experimental systems (White and Peterson 1990;White and Brantley 2003;Velbel 1993).
Another area which requires improvement is the development of strategies for coupling fluid flow and reactive transport. Global implicit, sequential-iterative or sequential noniterative schemes are often used to solve most reactive transport problems (Steefel et al. 2015). Irrespective of the numerical algorithm used, permeabilities and reactivities are updated with a time lag. Although this is not often done for porosity and diffusion coefficients, it does suggest that there is the possibility of errors in small mass balance at each time step. In addition, the simulation of an evolving porous media may be complicated when attempting to update the relevant parameters of the media as the evolution proceeds because of the interactions between flow and transport. In many published models, both multi-component and groundwater flow problems are solved sequentially. This does not apply, however, to evolving porous media because the evolution of porosity affects the flow equation through the storage term and modifies the permeability. Decoupling flow and reactive transport should help avoid numerically induced minor mass balance errors. Interestingly, there have been recent attempts to recouple both multiphase flow and reactive transport through the variation of porosity (Seigneur et al. 2018).
Processes occurring at the pore and macroscopic scales cannot be adequately described by porosity evolution alone. Unfortunately, the rigorous mathematical frameworks required for such an integrated approach are not yet available. Instead of these formulations, the geometrical organisation of the pore scale may be adequate. For instance, fluid-skeleton interactions could be described based on the surface area as pathways with high tortuosity. This is known to have lower transport properties caused by increased surface area.
Finally, In PHREEQC, the 1D advection-dispersion transport is calculated by an explicit finite difference algorithm. Sometimes, when the grid is coarse, the algorithm will compute numerical dispersion, which can be large or small in magnitude depending on the modelled reactions. Likewise, for the sorption reaction, the software has some uncertainties in defining the composition of the sorbed species, the number of sites, and the surface area. Therefore, most models of surface complexation applications depend on experimental data for specific materials from the studied site. PHREEQC also has a limited capability to integrate a system with complex kinetic reactions due to its explicit time simulation method (MacQuarrie and Mayer 2005).