A Review of Geochemical Modeling for the Performance Assessment of Radioactive Waste Disposal in a Subsurface System

: Radionuclides are inorganic substances, and the solubility of inorganic substances is a major factor affecting the disposal of radioactive waste and the release of concentrations of radionuclides. The degree of solubility determines whether a nuclide source migrates to the far ﬁeld of a radioactive waste disposal site. Therefore, the most effective method for retarding radionuclide migration is to reduce the radionuclide solubility in the aqueous geochemical environment of subsurface systems. In order to assess the performance of disposal facilities, thermodynamic data regarding nuclides in water–rock systems and minerals in geochemical environments are required; the results obtained from the analysis of these data can provide a strong scientiﬁc basis for maintaining safety performance to support nuclear waste management. The pH, Eh and time ranges in the environments of disposal sites cannot be controlled, in contrast to those under experimental conditions in laboratories. Using a hypothetical error mechanism for the safety assessment of disposal sites may engender incorrect assessment results. Studies have focused on radionuclide reactions in waste disposal, and have offered evidence suggesting that these reactions are mainly affected by the geochemical environment. However, studies have not examined the thermodynamics of chemical reactions or interactions between water and minerals, such as the surface complexation and adsorption of various nuclide-ion species. Simple coefﬁcient models have usually been applied in order to obtain empirical formulas for deriving Kd to describe nuclide distributions in the solid or liquid phase in water–rock geochemical systems. Accordingly, this study reviewed previous research on the applications of geochemical models, including studies on the development of geochemical models, sources of thermodynamic databases (TDBs) and their applications in programs, the determination of the adequacy of TDBs in surface complexation models and case studies, and the selection and application of activity coefﬁcient equations in geochemical models. In addition, the study conducted case studies and comparisons of the activity coefﬁcients derived by different geochemical models. Three activity coefﬁcient equations, namely the Davies, modiﬁed Debye–Hückel, and Pitzer equations, and four geochemical models, namely PHREEQC, MINEQL+, MINTEQA2, and EQ3/6, were used in the study. The results demonstrated that when the solution’s ionic strength was <0.5 m, the differences in the activity coefﬁcients between the Davies and modiﬁed Debye–Hückel equations were <5%. The difference between the Pitzer and Davies equations, or between the Pitzer and modiﬁed Debye– Hückel equations in terms of the calculated activity coefﬁcients was <8%. The effect of temperature on the activity coefﬁcient slightly inﬂuenced the modeling outputs of the Davies and modiﬁed Debye–Hückel equations. In the future, the probability distribution and uncertainty of parameters of Kd and the equilibrium constant can be used in geochemical and reactive transport models to simulate the long-term safety of nuclear waste disposal sites. The ﬁndings of this study can provide a strong scientiﬁc basis for conducting safety assessments of nuclear waste disposal repositories and developing environmental management or remediation schemes to control sites marred by near-surface contamination.


Introduction
A nuclear regulatory authority must provide a permit for the construction and operation of a repository for radioactive waste disposal. Therefore, a comprehensive multidisciplinary approach to safety assessment is required in order to determine whether the repository can safely dispose of radioactive nuclear waste. Currently, radioactive waste disposal often employs a multibarrier system composed of engineering barriers and natural barriers. The engineering barrier system includes the waste form, waste package, and concrete; furthermore, a buffer material is combined with backfill material to act as a physical and chemical barrier to retard the migration of radionuclides. The multibarrier system is designed for the isolation of radioactive waste from the biosphere. The evolution of repositories for geological disposal involves highly complicated coupling processes, including the mechanical stress change (M) of the repositories during excavation and construction, radionuclide decay heat (T) after the repository's closure, groundwater pressure (H) in the geological environment, and the interactions of various geochemical processes (C) among the groundwater chemical species, rock minerals and possibly released nuclides. The chemical field plays a crucial role after the disposal site's closure. Chemical changes occur because of near-field geochemical interactions caused by changes in temperature, the flow field and the stress field. Geochemical evolution affects the chemical composition of groundwater, waste form dissolution and nuclide release. These changes are considered nuclide migration problems in the performance safety assessment of repositories.
Therefore, geochemical models are required during the assessment process to predict the concentrations and chemical speciation of radionuclides entering the environment through different events and processes. A geochemical model is used to describe the water-gas-rock interaction in the aquatic system. The activity calculation method in a geo-chemical model is based on the ion-association theory. When the ionic strength of an aqueous solution increases to a value between 0.5 and 1 mol/kgw, the geochemical model can obtain reliable results. However, when the ionic strength exceeds this limit, a geochemical model based on an ion-interaction theory (such as the Pitzer equation) must be invoked for accurate predictions. There are two kinds of numerical methods for thermodynamic equilibrium modeling that can be applied to the calculation of speciation in aqueous systems, namely (1) the Law of Mass Action and (2) the Gibbs Energy Minimization. In general, two methods can be used to calculate the species distribution from a thermodynamic database. The first method consists in solving a set of nonlinear equations generated by the system's mass balances and equilibrium constants. The second method involves determining the most thermodynamically stable state by minimizing the reaction's free energy (the lowest energy state). Both methods are based on establishing the mass balance and chemical equilibrium state of the system. The relationship between the equilibrium constants, free energy, and the related chemical reactions and thermodynamic parameters is described below [1].
The following mass-action law defines an equilibrium constant in a thermodynamic database of a geochemical model: where a, b, c and d are the stoichiometric coefficients of the reactants A and B, and the products C and D, respectively, for the reaction (1). The terms a A , a B , a C and a D are the activities of A, B, C and D, respectively. K denotes the thermodynamic equilibrium constant.
In general, K can be characterized in the following way for these specific types of reactions: the complexation constant (K) for complex formation/dissolution reactions, the stability constant (K) for redox reactions, the solubility-product constant (Ksp) for dissolution/precipitation reactions, and the distribution coefficient (Kd) and selectivity coefficient (Kx) for sorption reactions. Surface complexation theory describes sorption phenomena based on the surface reaction equilibrium. Such phenomena occur on the surface of materials such as iron, aluminum, silica and manganese hydroxides, as well as on those of humic substances. Accordingly, four models have been proposed, namely the constant capacitance model, the diffuse double-layer model, the triple-layer model, and the charge distribution multisite complexation model. Surface complexation models are still being developed for the description of the adsorption/desorption reactions of nuclides. For example, Marmier [2] reported that cesium can be adsorbed on the surface of silica, and that a surface complexation model (SCM) can describe the adsorption of cesium on the surface of silica. Moreover, BarNes et al. [3] indicated that cesium tends to bind to silica particles in cement hydration minerals. Iwaida et al. [4] observed that cesium can be adsorbed on the silicate chain (C-S-H) of cement. In addition to cesium nuclides, Wang and Turner [5,6] reported that the adsorption of radionuclides is a major concern in water-mineral geochemical reactions with regard to the safety assessment of radioactive waste. They have developed a surface complexation model based on diffusion layer theory to assess the cations (such as Am 3+ , Pu 4+ , Pu 5+ and Np 5+ ) and anions (such as I − , IO 3 − , SeO 3 2− , SeO 4 2− and TcO 4 − ) of radionuclides in clay and various minerals by using chemical thermodynamic and equilibrium constant data.
For equilibrium reactions, the equilibrium constant is temperature-dependent, and can be calculated using the Van't Hoff equation, as follows: [1] K T = K Tr exp ∆H Tr R T − T r TT r (3) where K T and K Tr are the equilibrium constants for equilibrium reactions at the temperature of interest, T, and the referenced temperature, T r (298.15 K), respectively; ∆H Tr , is the enthalpy of the associated equilibrium reactions; and R is the ideal gas constant (8.315 J mol −1 K −1 ). The equilibrium constant is also pressure-dependent, and can be expressed as follows: lnK P = lnK P rS − ∆V T,P rS T · R · β ln ρ P ρ P rS (4) where K P and K P rS are the equilibrium constants for the equilibrium reactions at the pressure of interest, P, and the reference pressure, P rs (i.e., the saturated vapor pressure), respectively; ∆V T,P rS is the volume change of the dissociation reaction at the temperature T and reference pressure P rs ; β is the isothermal compressibility coefficient of water at the temperature T and the pressure of interest, P; and ρ P and ρ P rS are the densities of water at the pressure of interest, P, and reference pressure, P rs , respectively. [1] The equilibrium constant can also be related to the standard molal Gibbs free energy of a reaction at any temperature and pressure by the following equation: ∆G 0 r = ∆G 0 r,P r ,T r −∆S 0 r,P r ,T r (T − T r )+ T T r ∆C 0 r,P r dT − T T T r ∆C 0 r,P r dlnT+ P P r ∆V 0 r dP (6) where T r and P r stand for the reference temperature and pressure, and ∆G 0 r,Pr,Tr , ∆S 0 r,Pr,Tr , ∆C 0 r,Pr,Tr and ∆V 0 r are the stand molal Gibbs free energy, entropy, heat capacity and volume of the reaction at the temperature and pressure indicated by the subscripts, respectively [7].
In the reaction (1), the quantities of the reactants A and B, and the products C and D are represented as activities, a i , and not as concentrations, C i , with respect to a species, i. The activity of a species is related to its concentration at a specific ionic strength by the following equation: The activity coefficient, γ i , is an ion-specific correction factor describing the interactions among charged ions in an aqueous solution. Because the activity coefficient is a nonlinear function of ionic strength given by the Davies, modified Debye-Hückel, or Pitzer equations, the activity is also a nonlinear function of concentration.
When the ionic strength reaches 0.1 mol/kg, the activity of a chemical species in aqueous solution decreases due to the interaction among the charged ions. The activity decreases with an increase in the ionic strength, and is always lower than the concentration. Therefore, the activity coefficient is less than 1. It can be concluded that the decrease in activity will be more pronounced at higher ion concentrations and higher ion valences. In an infinitely dilute solution, the interaction between ions will approach zero; therefore, the activity coefficient will be 1, and the activity will be equal to the concentration [1].
In Raoult's law, the conditions of an ideal solution are based on the assumption that chemical interactions between bonded molecules are equal, but the interactive forces between ions are zero. However, in real solutions, attractive forces such as adhesive and cohesive forces are not uniform between the molecules, thus causing real solutions to deviate from Raoult's law. Consider the effect of ion interaction on the activity of an electrolytic solution; the ionic strength, I, represents the collective concentration level of all of the ions in the solution.
where C i is the molal concentration of ion i (mol/kgw), Z i is the charge number of ion i, and Σ is summation of the charge numbers of all of the ions in the solution.
In general, an activity coefficient varies with the concentration of a solution. Different activity coefficient models may yield significantly varied activity coefficients. In this study, we first reviewed the theory of the activity coefficient and compare the applicable concentration range of the activity coefficient models.
1. Debye-Hückel Theory The Debye-Hückel theory [8] states that the activity of a solution is affected by the mutual electrostatic forces between the ions. Thus, the activity coefficient γ is introduced in order to calculate the correct activity of an ion in a solution, and the Debye-Hückel extended equation can be expressed as follows: where A and B are temperature-dependent constants, and a i is the diameter of ion i. The Debye-Hückel limiting law is simplified to (3) in a dilute solution.
The Debye-Hückel equation is commonly used to estimate the activity coefficient of a dilute single-ion solution. The Davies equation and modified Debye-Hückel equation are used to describe complicated ionic interactions in high-concentration solutions, and these equations are extensively applied to determine activity coefficients at different concentrations.

of 35
The Davies equation [9] is an empirical extension of the Debye-Hückel equation, which can be used to calculate the activity coefficients of electrolyte solutions at a relatively high concentration, as follows: where b is an empirical parameter ranging from 0.1 to 0.4. In the EQ3/6 model, the constant b is assigned a fixed value of 0.2, the PHREEQC and MINTEQA2 are 0.3, and MINEQL+ is set to 0.24 [10]. The modified Debye-Hückel equation (D-H), also called the B-dot equation [11,12], is an improvement over the Debye-Hückel equation, and is presented in (12).
where * B is a temperature-dependent ion-specific parameter, and a i is the theoretical distance of the closest approach between ions and opposite charges, but in practice, it is an adjustable parameter.
Considering the effect of short-range interactions between ions and hydration reactions, the Davies model and modified Debye-Hückel (B-dot) model derived using polynomials of ionic strength can predict activity coefficients more precisely than the Debye-Hückel model when the ionic strength increases (Figures 1 and 2).    Comparisons of the activity coefficients for Ca 2+ and Cl − with ionic strength (CaCl 2 solution, a i = 4.86 for Ca 2+ , a i = 3.71 for Cl − [13]) are shown in Figures 1 and 2, respectively. The figures demonstrate the substantial differences in the activity coefficients calculated for calcium and chloride based on the Debye-Hückel, Davies and B-dot models. The predictions clearly show that different ions have different activity coefficients at the same ionic strength.
3. Pitzer equation The Pitzer Equation [14], developed for deriving the excess Gibbs free energy of electrolytic solutions, combines the Debye-Hückel equation with additional terms in the form of a virial equation, and it is widely used in geochemistry to model systems of mixed or single strong electrolytes at high concentrations, such as seawater.
where γ ± is the mean activity coefficient; ν M and ν X are the numbers of M and X ions in the formula, respectively; Z M and Z X are their respective charges in electronic units; and m is the molality. The other quantities are defined as follows: where A φ is the Debye-Hückel coefficient for the osmotic function, B γ MX is the second virial coefficient, and C γ MX is the third virial coefficient. When a radioactive waste repository is in a brackish water area, the surrounding groundwater can have a high ionic strength. If the concentration of a high-salinity ground-Appl. Sci. 2021, 11, 5879 7 of 35 water quality is obtained by laboratory analysis, it can be difficult to determine the real reaction concentration values of the groundwater quality due to the chemical interactions among all of the relevant ion species. The interaction among the ionic species in groundwater is related to the concentrations of the cations and anions in the aqueous solution, and the charge of each species. Therefore, it is necessary to analyze the influence of the ionic strength and the activity coefficients on the activity of the chemical species. These calculations will allow for the accurate prediction of the speciation and solubility of the nuclides released from the repository to the groundwater, and the sorption capacity of the nuclides adsorbed by bentonite.
Geochemical simulation can serve as an essential tool in the support of repository performance for two main reasons. First, geochemical simulation can help identify the chemical species formed and predict the future geochemical/environmental state. Second, there can be different scenarios and geochemical processes affecting the safety of a repository. Numerical geochemical modeling offers an effective (and often the only viable) means to conduct simulations assessing the performance of radioactive waste disposal for differing scenarios and conditions. With this purpose in mind, this study reviewed the literature on geochemical models and their applications. The studies include geochemical models, the sources of thermodynamic databases (TDBs) used in models, using TDBs to simulate nuclide solubility, the 'sufficiency' of TDBs for surface complexation models, case studies of nuclide sorption, and the selection and application of activity coefficient equations in geochemical models. In addition, case studies were conducted, and the performance of different activity coefficient equations in the geochemical models was compared.

Development of Geochemical Models
Geochemical models enable the conceptualization of chemical reactions to be represented by mathematical equations, which can then be coded up as a geochemical modeling software tool. The first geochemical model was introduced by Garrels and Thompson in 1962 [15]. Garrels and Thompson established a chemical model of seawater based on the major ions present in seawater. In the 1980s, researchers studied water-rock interactions from the perspective of physical chemistry and further integrated solute transport, dispersion transport, chemical reactions and water-rock interface chemistry. A coupled model of solute transport was therefore developed. Over the last few decades, geochemical models have been developed rapidly in response to an increasing number of problems in geochemical and environmental science which require quantitative analysis and assessment. In research supporting nuclear waste disposal, geochemical models have become a crucial tool for the simulation of geochemical processes, including the formation and evolution of hydrochemical components in water-rock systems, and surface complexation between minerals and nuclides. Geochemical models can simulate many geochemical phenomena and processes, covering: (1) chemical species in water; (2) the saturation index (SI), which can be expressed as SI = log(IAP/Ksp), where IAP is the ionic activity product, and Ksp is the solubility product; (3) the mixing of more than two solutions, as well as the resulting aqueous species composition, Eh, pH, and SI; (4) dissolution (precipitation) caused by the interaction of organic and inorganic compounds; (5) reduced pressure or changes in water's chemical composition and SI caused by boiling or the addition of a gas; (6) the rates or kinetics of geochemical reactions; and (7) water-rock interactions. Simulations of water-rock interactions can be further divided into forward and inverse modeling. Forward modeling involves the simulation of possible dissolution-precipitation, adsorption-desorption and oxidation-reduction reactions, and the prediction of the chemical species involved in the dissolution or precipitation, as well as the mass being adsorbed or desorbed by the minerals. Inverse modeling involves the simulation of the aqueous species composition and Eh and pH values in water-rock interactions.
As chemical effects are highly crucial in the closure of a high-level radioactive waste disposal site, a geochemical model can provide an effective means to calculate the formation and evolution of chemical species under long-term water-rock interactions in both the engineering and natural barriers of the repository. Such a model is usually based on the thermodynamic data of chemical reactions. The geochemical reaction is represented by mathematical equations, which can be written into a numerical model and solved using a mathematical algorithm for its compilation into computer code. The approach of the law of mass action requires the concentration of the master chemical component and thermodynamic data such as the equilibrium constant of the product species at the pressure and temperature of interest in the heterogeneous system. For the Gibbs energy minimization, a value of the standard Gibbs free energy of each chemical component is taken from the thermodynamic database and corrected to the pressure and temperature of interest, if necessary.
Amongst the most widely used geochemical models in the approach of the law of mass action are SOLINEQ.88 [16], EQ3/6 [10,[17][18][19], MINEQL+ [20], MINTEQA2 [21], PHREEQC [22] and WATEQ4F [23]. The represented geochemical models for the approach of Gibbs energy minimization are GEM Selektor [24], Reaktoro [25,26], ChemApp [27] and ChemSage/FactSage [28,29]. Geochemical models can be used for the safety assessment of engineering and natural barriers in high-level radioactive geological disposal repositories, and they can serve as a research tool for the determination of water-rock interactions and the reactive chemical transport of nuclides in near and far fields. These models are mainly composed of several modules, such as input, output, iterative algorithm and chemical TDB modules. In general, the input module consists of the total concentration of the components in the initial solution as well as the solid and gas components. The output module primarily includes the equilibrium components, species concentrations and mass exchange between the aqueous solution, solid phase and gas phase. The iterative algorithm module is used to solve nonlinear equations, and it often employs the improved Newton Raphson iterative method. The database module of chemical thermodynamics includes data on the thermodynamic equilibrium constant under a standard state (25 • C, 1 atm (1 × 10 5 Pa)), the equilibrium constant at different temperatures, the parameters required for activity coefficient calculation, and the stoichiometric factors of each species in aqueous solutions.

Evaluation of a Thermodynamic Database (TDB)
The geochemical models are shipped with thermodynamic databases, and each equilibrium constant of the TDBs for a chemical reaction can be separated from the simulation model. This means that the equilibrium constant is not rigidly embedded in the simulation model. Users can easily change the equilibrium constant value and add/delete chemical reactions without affecting the code function in the database. The geochemical models used in this study (such as PHREEQC and EQ3/6) allow users to modify the equilibrium constants in input files or databases with internal consistency.
The required data of the equilibrium constant or standard Gibbs free energy can be found in the thermodynamics database of the geochemical model. We should check whether the database contains the thermodynamic and kinetic parameters of the species involved in the model. At the same time, we should focus on whether the data in the thermodynamic database have internal consistency; for example, the enthalpy of formation (∆ f H 0 ), standard Gibbs energy of formation (∆ f G 0 ), and entropy (∆ f S 0 ) for the chemical species x may come from different sources. The value (i.e., equilibrium constant or standard Gibbs free energy) of thermodynamics database may be obtained through the statistical regression analysis of experimental data at the beginning of the derivation; however, the value may not be consistent with the relationship of the formula (∆ f G 0 . Achieving the consistency of these parameters is an important task that is still in progress in the establishment of geochemical models and chemical thermodynamics databases [30,31]. The widely used PHREEQC geochemical model already has an internal TDB. In addition to the databases in the PHREEQC, WATEQ4F and MINTEQ models, Lawrence Livermore National Laboratory in the United States developed a database of chemical ther-modynamics. The database comes with EQ3/6, which can be used in the PHREEQC model. Andra/RWM, France, developed a chemical TDB named ThermoChimie-TDB [32]. To date, the format of the ThermoChimie database has been employed by other models, including PHREEQC, CHESS [33], Crunch [34] and Toughreact [35]. The compatibility of the TDBs can be extended to other models, including Medusa [36] and Geochemist's workbench [37], according to the demand for databases. Other research institutions have also developed various chemical TDBs. Examples include the Thermochemical Database Project developed by the Nuclear Energy Agency (NEA) [38][39][40]; the Thermodynamic Reference Database (THEREDA) focusing on transuranic elements, which was developed by various German institutions [41]; Japanese JNC TDB [42]; British HATCHES [43]; Nagra/PSI chemical TDB, developed by the Swiss Radioactive Waste Management Agency [44]; and MOLDATA [45] TDB, which is under compilation by SCK.CEN at Mol (Belgium) for use in a framework of safety case development by ONDRAF/NIRAS. These TDBs have been developed to support the geochemical modeling required for the safety and performance assessment of radioactive waste repositories. Together, they provide good coverage of geochemical reaction data for the hydrogeochemical simulation of radioactive waste disposal. Three TDBs-ThermoChimie-TDB, THEREDA and NEA TDB-are described as follows.
The latest release of ThermoChimie (V10a) can be downloaded using the "Download" table. The database contains the basic major elements, numerous nuclides and toxic chemical elements. The main elements include C, N, O, F, Na, Mg, Al, Si, P, S, Cl, K, CA, Fe and Ba. The nuclides and toxic chemical elements are outlined as follows: (1) long half-life nuclides and actinide daughter nuclides over 1000 years, including C, Al, Cl, K, Ca, Mn, Fe, Ni, Se, Rb, Zr, Nb, Mo, Tc, Pd, Sn, I, Cs, Ba, Sm, Ho, Hf, Pb, Ra, Th, Pa, U, Np, Pu, Am and Cm (Be, La, Bi and Ac will be included in the future); (2) short half-life nuclides, including CO, Sr, Ag and Eu; and (3) toxic chemical elements, including B, Cr, As, Cd and Sb. The basis for the development of the ThermoChimie database includes a pH range of 5-14, the redox stability range of water, a temperature of 15-90 • C, and an ionic strength between very low (0.001 mol/kgw) and 3 mol/kgw. The data selection process must consider both temperature and ionic strength correction. In many cases, there is not enough data to perform these corrections, and strategies have been developed to overcome this limitation [32,46].
The purpose of the THEREDA chemical database is to simulate the geochemical evolution of radioactive and chemical toxic waste disposal sites. In the process of modeling, it is necessary to check the internal consistency of the thermodynamic data (including the adsorption data) of chemical reactions between nuclides and the surrounding rock. The built-in databases provided by different geochemical models are inconsistent, such as the activity coefficient equation used in the data generation process and the experimental data of the equilibrium constant. Therefore, using different chemical databases to compare geochemical simulation results may have different results. This is due to the incomplete and conflicting data in the existing chemical databases, the different validities of the solution composition (ionic strength) and model limitations, and the lack of sorption data. THEREDA is a collaborative project that has addressed these challenges. The aim of the project is to establish a consistent and quality-assured database for all of the safetyrelevant nuclide and toxic element, temperature and pressure ranges, with a focus on high-concentration brine systems. The database uses the Pitzer equation to calculate the activity coefficient applicable to such conditions. The chemical elements and nuclides in the THEREDA database are as follows: H, C, O, Na, Mg, Al, Si, P, S, Cl, K, Ca, Se, Sr, Tc, Cs, Nd, Th, U, Np, Pu, Am and Cm. Part of the THEREDA project involves the collection and evaluation of uranium chemical reaction data. The aquatic chemical and thermodynamic conditions for U (VI) and U (IV) are crucial for the geochemical modeling of disposalsite-relevant systems. The THEREDA chemical database established the correction of the activity coefficients for U, Na, K, Mg, CA, Cl, S, and C in high-salinity aqueous solutions using the Pitzer equation. The THEREDA database does not have the same strict quality requirements for thermodynamic data as the NEA TDB. Currently, THEREDA is jointly planned, implemented and further developed by the research institute engaged in the final disposal of radioactive waste and its potential host rock [41].
The NEA TDB has provided high-quality and accurate chemical thermodynamic data used in the safety assessment of deep geological repositories for the international radioactive waste management community. The scenario analysis and assessment of radionuclide migration paths to the environment are essential for the performance assessment of geological disposal. The safety assessment of radioactive waste disposal is based on the use of a chemical thermodynamics database for the numerical simulation of the speciation. The NEA mentioned [47] that their TDB has become the worldwide standard for geochemical modeling in the field of geological disposal. However, the applicability of the NEA TDB data is limited to low-to-medium ionic strength environments and temperatures of less than 100 • C. Additionally, the NEA TDB data do not cover all of the species and thermodynamic quantities required for safety assessment. Therefore, for application to natural environmental systems, users may need to obtain data from other thermodynamic databases; however, it is necessary to ensure that the selected end states are consistent and represent the chemical reactions of the system of interest. The chemical elements and nuclides in the NEA TDB are as follows: H, B, C, N, O, F, Na, Mg, Al, Si, P, S, Cl, K, Ca, Cr, Mn, Fe, Co, Ni, Cu, Zn, Ga, Ge, As, Se, Br, Rb, Sr, Zr, Tc, Ag, Cd, In, Sn, Te, I, Cs, Ba, Hg, Tl, Pb, Th, U, Np, Pu and Am. The NEA TDB has the following characteristics: a comprehensive literature review and detailed discussion of the performed evaluation; the traceability of all used the documents used in the evaluation and selection of the data; strict guidelines and a strengthened selection of high-quality values; and thermodynamic values conform to the CODATA standards [31,48,49].

Radionuclide Speciation and Solubility
The radionuclides selected for the geochemical simulation were Thorium-232, Neptunium-237 and Plutonium-242, which have a long half-life and are safety-relevant materials for radioactive waste in Taiwan. We referred to the test calculations from THEREDA [50] and used the thermodynamic database (THEREDA_PIT_PHRC_R-07.dat) to simulate the radionuclide speciation and solubility for the solid compounds Th(OH) 4 , Np(OH) 4 and Pu(OH) 4 in NaCl, CaCl 2 , Na 2 CO 3 and KHCO 3 aqueous solution conditions (Tables 1-3) [50] using the PHREEQC model. Additionally, we also used the PHREEQC model, and used the NEA TDB (NEA_TDB_phreeqc_Nov2018.dat) to simulate radionuclide speciation and solubility for the solid compounds Th(OH) 4 , Np(OH) 4 and Pu(OH) 4 in the groundwater (Table 4) of subsurface systems of potential radioactive waste disposal sites in Taiwan.
The simulation results showed that the main dissolved ions of Th(OH) 4(s) are Th 4+ (aq) and Th(OH) 4(aq) , and the total dissolved concentration of Th 4+ can vary with the pH, as shown in Figure 3. and Pu 4+ (aq) . The total dissolved concentration of Pu 4+ also varied with the pH ( Figure 5). There are 16 aqueous complexation formations for Th(OH) 4(s) , 30 aqueous complexation formations for Np(OH) 4(s) , and 24 aqueous complexation formations for Pu(OH) 4(s) in the dissolution reaction of a solid compound. Figure 6 showed the total dissolved concentrations of the solid compounds Th(OH) 4 , Np(OH) 4 and Pu(OH) 4 , as calculated using the groundwater quality of a potential radioactive waste disposal site in Taiwan, which can be used as the initial concentration data of radionuclides in the adsorption/desorption interactions of surface complexation reactions between nuclides and bentonite after its release from the waste canister to the engineer barrier.  The simulation results showed that the main dissolved ions of Th(OH)4(s) are Th 4+ (aq) and Th(OH)4(aq), and the total dissolved concentration of Th 4+ can vary with the pH, as shown in Figure 3. The main dissolved ions of NpCl 3+ (aq), Np(OH) 3+ (aq), Np(OH)2 2+ (aq), Np(OH)3 + (aq), Np(OH)4(aq) and Np 4+ (aq) have a valence of 4, and the redox formations NpO2 + (aq), NpO2(OH)(aq) and NpO2(OH)2 − (aq) have a valence of 5. The total dissolved concentration of Np 4+ in relation to the pH is shown in Figure 4. The main dissolved ions of Pu(OH)4(s) are PuCl 3+ (aq), Pu(OH)4(aq), Pu(OH)3 + (aq), Pu(OH)2 2+ (aq), Pu(OH) 3+ (aq) and Pu 4+ (aq). The total dissolved concentration of Pu 4+ also varied with the pH ( Figure 5). There are 16 aqueous complexation formations for Th(OH)4(s), 30 aqueous complexation formations for Np(OH)4(s), and 24 aqueous complexation formations for Pu(OH)4(s) in the dissolution reaction of a solid compound. Figure 6 showed the total dissolved concentrations of the solid compounds Th(OH)4, Np(OH)4 and Pu(OH)4, as calculated using the groundwater quality of a potential radioactive waste disposal site in Taiwan, which can be used as the initial concentration data of radionuclides in the adsorption/desorption interactions of surface complexation reactions between nuclides and bentonite after its release from the waste canister to the engineer barrier.

Evaluation of the Thermodynamic Data for the Surface Complexation
The thermodynamic data of SCM can be used by geochemical and reactive transport models to evaluate and compare the composition and concentration of radionuclides in water and minerals in geochemical reactions in various hydrogeochemical environments. Bradbury et al. [51] proposed a surface complexation model to simulate radionuclide adsorption in natural systems. They indicated that the results of experimental methods, such as distribution coefficients (Kd) and isothermal equations, are valid only under experimental measurement conditions, and that the application of such results to other conditions and systems may be problematic. Additionally, they studied the diffuse double-layer model to describe the adsorption of radionuclides by igneous minerals. They reported that surface complexation adsorption, the development of a conceptual model, and the use of a prediction model constitute a feasible method for the scenario analysis of the safety assessment of geochemical reactions in radioactive waste environments. This method can be used under the environmental conditions of disposal sites (varying pH, Eh, and ionic composition of the surrounding water). Surface complexation models can

Evaluation of the Thermodynamic Data for the Surface Complexation
The thermodynamic data of SCM can be used by geochemical and reactive transport models to evaluate and compare the composition and concentration of radionuclides in water and minerals in geochemical reactions in various hydrogeochemical environments. Bradbury et al. [51] proposed a surface complexation model to simulate radionuclide adsorption in natural systems. They indicated that the results of experimental methods, such as distribution coefficients (Kd) and isothermal equations, are valid only under experimental measurement conditions, and that the application of such results to other conditions and systems may be problematic. Additionally, they studied the diffuse doublelayer model to describe the adsorption of radionuclides by igneous minerals. They reported that surface complexation adsorption, the development of a conceptual model, and the use of a prediction model constitute a feasible method for the scenario analysis of the safety assessment of geochemical reactions in radioactive waste environments. This method can be used under the environmental conditions of disposal sites (varying pH, Eh, and ionic composition of the surrounding water). Surface complexation models can provide reliable results for safety assessments in radioactive waste disposal.
The Nuclear Energy Agency of the Organization for Economic Co-operation and Development has regularly discussed and studied geochemical reactions concerning the surface complexation adsorption of nuclides. The agency has executed such discussions and studies in collaboration with research organizations such as ANDRA and IPSN, France; BNFL and NIREX, the United Kingdom; ENRESA, Spain; KAERI, Korea; NAGRA, Switzerland; ONDRAF/NIRAS, Belgium; JNC, Japan; POSIVA, Finland; PSI, Switzerland; SKB and SKI, Sweden; and RAWRA, Czech Republic. The Kd value of such reactions is estimated and used in the safety and performance assessment of radioactive waste disposal sites [52].

Application of Thermodynamic Data in SCM
The total dissolved concentration of Np 4+ (Figure 2) was used as the concentration of the nuclide released from the waste disposal site. PHREEQC was used to simulate the surface complexation reaction of the neptunium-237 radionuclide on montmorillonite for the groundwater quality (Table 4) in the subsurface system of the potential disposal site in Taiwan. MX-80 bentonite is a buffer material used in radioactive waste disposal facilities; it is composed of 65-75% montmorillonite. Bradbury et al. [53] developed equilibrium constants for the surface complexation reactions of the sorption of Np 4+ nuclides on MX-80 bentonite ( Table 5). The equilibrium constants of these surface complexation reactions were used as the parameters of the chemical reactions in this study. No currently available database contains the SCM parameters that are defined in the NEA TDB (NEA_TDB_ phreeqc_Nov2018.dat) that we used in the PHREEQC modeling. Therefore, the thermodynamic data for Np nuclide adsorption to montmorillonite (Table 5) were included in the input file of the PHREEQC. With a complete database of chemical thermodynamics, we can use the PHREEQC to simulate the complexation reaction of Np 4+ nuclide sorption on montmorillonite and other geochemical reactions for the groundwater quality of the potential disposal site in Taiwan.
The results showed that neptunium with a valence of 3 forms various chemical species due to redox reactions. However, there was no solid nuclide precipitation, only nuclide adsorption on montmorillonite. Neptunium with a valence of 3 forms trivalent, tetravalent, pentavalent and hexavalent ions. The main dissolved ions are as follows: trivalent neptunium, i.e., Np(CO 3 ) 3 3− , NpOH 2+ , and Np 3+ ; tetravalent neptunium, i.e., Np(OH) 4 6 6− . The results also showed that the value of Kd is 36.7 cm 3 /g at a pH of 8.996 and a pe of 4.068 according to the adsorption capacity of the surface complexation reaction and the total dissolved concentration ( Table 6).

Different Activity Coefficient Equations in Geochemical Models
Following groundwater invasion into a deep geological repository, radioactive species are released and react with hydrogeological materials. Considering the ion interactions occurring in the in-situ chemical reactions, an activity coefficient is used to modify the chemical activities to represent the real concentration of the generated species.
Four chemical equilibrium models, namely PHREEQC, MINEQL+, MINTEQA2 and EQ3/6, were used to test the model's applicability. The results can provide a general guideline for the selection of an appropriate model in field assessments.

Geochemical Equilibrium Models
Modeling the chemical equilibrium of target compounds is relevant during the control of pH, alkalinity, temperature and other factors. As part of this approach, one must determine the rates of the contribution of various components with specific chemical characteristics to mixtures. A model established using thermodynamic equilibrium theory and experimental data is suitable for the probing of the reaction of geochemical equilibrium; furthermore, the model can present the reaction process in a natural environment in numerical form. With extensive theoretical research and experimental data analyses on geochemical reactions, geochemical equilibrium models have become increasingly reliable for the simulation of in situ subsurface water-rock interactions and the prediction of solution concentrations through mineral ingredient analysis. Thermodynamic chemical equilibrium typically aims at the determination of spontaneous changes to an equilibrium state. Any system that is not at equilibrium would change spontaneously with the release of energy. In general, simulation is more straightforward for a closed system than an open system, in which the equilibrium condition is simplified with mass balance and thermodynamic equilibrium. This study compared the performance of four chemical equilibrium models, namely PHREEQC, MINEQL+, MINTEQA2, and EQ3/6, in terms of the prediction of activity coefficients.

PHREEQC Model
PHREEQCI is a complete windows-based graphical user interface to the geochemical computer program PHREEQC developed by the U.S. Geological Survey, and it provides all of the capabilities of the geochemical model PHREEQC, including speciation, batchreaction, 1D reactive-transport and inverse modeling [22].

MINEQL+ Model
MINEQL+ is a chemical equilibrium model capable of calculating aqueous speciation, solid phase saturation states, precipitation-dissolution and adsorption. An extensive thermodynamic database is included in the model [20].

MINEQA2 Model
MINEQA2 was originally developed by Battelle Pacific Northwest Laboratory in the mid-1980s as MINTEQ. MINTEQA2 was derived from MINTEQ. MINTEQA2 is a chemical equilibrium model for dilute solutions with ionic strengths up to about 0.5 m. The model is useful for the calculation of the equilibrium mass distribution among dissolved species, adsorbed species and multiple solid phases under a variety of conditions, including a gas phase with constant partial pressures. A comprehensive database is included, which is adequate for solving a broad range of problems without the need for additional usersupplied equilibrium constants [21].

EQ3/6 Model
In the mid-1970s, the EQ3/6 software package was developed at Northwestern University to model seawater-basalt interactions in mid-ocean ridge hydrothermal systems. The EQ3/6 contains two kinds of calculations for aqueous solutions and aqueous systems. The first kind is called a speciation-solubility calculation, which describes the chemical and thermodynamic state of the solution using input analytical data and theoretical assumptions. The second kind of calculation is called a reaction path calculation. A reaction path calculation predicts the evolution of a reacting system [18,19].
The four chemical equilibrium models provided several equations relevant to activity coefficient calculation (Table 7). Table 7. Options of activity coefficient equations incorporated with four models.

Model Calibration
According to the analysis by Nordstorm et al. [54], the simulation of different chemical models in the same conditions brings about different outcomes due to the different thermodynamic databases in the models. The thermodynamic database includes the equilibrium constant, heat constant and stoichiometry. In the database of the geochemical model, there are chemical species that are mathematically independent, which are generally called chemical components. The chemical components in the database of the geochemical model are the main cation and anion species of in-situ water quality data. We can select these chemical components as reactants in order to calibrate the model. Thus, in order to realize the effect of the database on the chemical equilibrium model, a test case with the same concentrations of all of the reactants from the in-situ water quality data was selected as the input of the four models (Table 8). The concentration and activity coefficients were calculated using the Davies equation with their respective thermodynamic databases ( Table 9).
The concentration of the reaction products (Table III) and activity coefficients (Table IV) given by the four models demonstrate great differences. Figures 7 and 8 show the concentration and activity coefficient of the major reaction products calculated by the Davies equation using the built-in thermodynamic databases of the four models, respectively. The concentrations of the reaction products and activity coefficients obtained from the four models demonstrate great inter-model differences, especially for MINEQL+. The reason for this is that MINEQL+ produces precipitation species, such as dolomite and calcite, while the saturation index (SI) values calculated by the other three model systems show that precipitation reactions do not easily occur. In order to remove the difference caused by these databases on the activity coefficients, we took the "phreeqc.dat" database of PHREEQC as the standard in order to revise the databases MINTEQA2 and COM, and recalculated new concentration of products ( Figure 9) and activity coefficients ( Figure 10). The results showed that only the concentrations of Na + , NaHCO 3 and NaCO 3 − in the EQ3/6 model were different from those in the other models. The reason for these differences is that there are chemical reactions of NaCl and NaOH in the COM database. These reactions were not considered in the PHREEQC and MINTEQA2 databases. Thus, the results showed that the differences of the concentration and activity coefficient between the models decreased significantly (Figures 7-10) and assured us that the use of the same thermodynamic database can produce similar activity coefficients through the four models.  the "phreeqc.dat" database of PHREEQC as the standard in order to revise the databases MINTEQA2 and COM, and recalculated new concentration of products ( Figure 9) and activity coefficients ( Figure 10). The results showed that only the concentrations of Na + , NaHCO3 and NaCO3 -in the EQ3/6 model were different from those in the other models. The reason for these differences is that there are chemical reactions of NaCl and NaOH in the COM database. These reactions were not considered in the PHREEQC and MINTEQA2 databases. Thus, the results showed that the differences of the concentration and activity coefficient between the models decreased significantly (Figures 7-10) and assured us that the use of the same thermodynamic database can produce similar activity coefficients through the four models.     At present, there is no thermodynamic database in any geochemical model that contains all of the chemical species and reaction types of the solutions and minerals that may react. Developers cannot guarantee the internal consistency of the TDB. Recent research has focused on the determination of the optimal TDBs and the internal consistency of the data (e.g., standard Gibbs free energy) [30,31]. Thus, the main challenges and limitations of geochemical modeling would be to check the accuracy of the data in the database, the consistency of the chemical thermodynamics data, and the thermodynamic data sensitivity of the simulation using the database.  At present, there is no thermodynamic database in any geochemical model that contains all of the chemical species and reaction types of the solutions and minerals that may react. Developers cannot guarantee the internal consistency of the TDB. Recent research has focused on the determination of the optimal TDBs and the internal consistency of the data (e.g., standard Gibbs free energy) [30,31]. Thus, the main challenges and limitations 0.

Comparison of the Activity Coefficients between the Models
This study considered various concentration levels in a solution (Ca 2+ , HCO 3− , Cl − , K + , Mg 2+ , Na + , SO 4 2− , Fe 3+ and PO 4 3− were considered to be the major ions) in order to investigate the differences in the activity coefficient calculation performance between the four models (i.e., PHREEQC, MINEQL+, MINTEQA2 and EQ3/6).
The considered cases were as follows: (1) The activity coefficient of a species calculated using the Davies equation in the four models.
(2) The activity coefficient of a species calculated using the modified Debye-Hückel equation in EQ3/6 and PHREEQC.
(3) The activity coefficient of a species calculated using the Pitzer equation in EQ3/6. (4) Considering the effect of temperature on the activity coefficient, the activity coefficient of a species calculated using the modified Debye-Hückel and Davies equations at 10 • C, 25 • C, 60 • C and 90 • C.

Calculation Results
The results showed that if the ionic strength (I) of a solution was <0.0001 m, the activity coefficients of any species corrected using the Davies or modified Debye-Hückel equations for different charge numbers of ions (Z) were mostly greater than 0.95 (Table 10). It should be noted that except for Fe 3+ ions with Z = 3, which are slightly less than 0.95, all other ions meet this threshold value (I < 0.0001 mol/kgw). The results suggest that the activity coefficient correction is small and can be neglected at such a low concentration level. When I was <0.5 mol/kgw, the activity coefficient variation calculated using the Davies equation was <3%; as I increased from 0.5 to 0.9 mol/kgw, the maximum variation in the activity coefficient increased to about 10% ( Figure 11). Therefore, the Davies equation is suitable for solutions with an ionic strength of <0.5 mol/kgw. Moreover, the activity coefficients of high-charge species calculated using the Davies or modified Debye-Hückel equations were more sensitive than those of low-charge species. Comparing the activity coefficients calculated using the Davies and modified Debye-Hückel equations through PHREEQC and EQ3/6 revealed that if I was <0.5 mol/kgw, the activity coefficient differences were <5%; if I was >0.5 mol/kgw, the activity coefficient differences increased slowly ( Figure 12). According to the theory of activity coefficient, the Davies equation can be suitable for solutions with I < 0.5 mol/kgw, and the modified Debye-Hückel equation can be suitable for solutions with I between 0.3 and 1 mol/kgw. If I was <0.5 mol/kgw, the differences in activity coefficients calculated using the Pitzer, Davies and modified Debye-Hückel equations through EQ3/6 were <8%. When I was >0.5 mol/kgw, the differences increased gradually ( Figure 13). Notably, the Pitzer equation described the activity coefficient satisfactorily for solutions with high ionic strength ( Figure 14). As the temperature increased from 10 • C to 90 • C, the activity coefficients calculated by the PHREEQC model using the Davies and modified Debye-Hückel equations decreased. Moreover, if I was <0.5 mol/kgw, the activity coefficients calculated using the Davies and modified Debye-Hückel equations were similar (Figure 15), suggesting that the Davies equation is applicable for conditions involving moderate temperature variations. The activity coefficient decreased with the increase of temperature from 10 • C to 90 • C for the modified Debye-Hückel equation, as calculated from EQ3/6 ( Figure 16).

Field Applications of the Two Sites
This study also considered two field cases, namely Shiau-Chiou in Kinmen County and Bai-Shi Lake in Hsinchu County. The four models were used to calculate the activity coefficients of the solutions at these two sites.
Shiau-Chiou is a possible disposal site for low-level radioactive waste in Kinmen County. The groundwater quality data measured by the Industrial Technology Research Institute (Table 11) were used in this study; the data were entered into the four models in order to calculate the activity coefficients of the solution ions using the Davies and modified Debye-Hückel equations. The ionic strength calculated for the groundwater samples obtained from Shiau-Chiou ranged from 0.3 to 0.5 mol/kgw (Figure 17), indicating that the Davies and modified Debye-Hückel equations can yield fairly consistent activity coefficients ( Figure 18). ature increased from 10 °C to 90 °C, the activity coefficients calculated by the PHREEQC model using the Davies and modified Debye-Hückel equations decreased. Moreover, if I was <0.5 mol/kgw, the activity coefficients calculated using the Davies and modified Debye-Hückel equations were similar (Figure 15), suggesting that the Davies equation is applicable for conditions involving moderate temperature variations. The activity coefficient decreased with the increase of temperature from 10 °C to 90 °C for the modified Debye-Hückel equation, as calculated from EQ3/6 ( Figure 16).

Field Applications of the Two Sites
This study also considered two field cases, namely Shiau-Chiou in Kinmen County and Bai-Shi Lake in Hsinchu County. The four models were used to calculate the activity coefficients of the solutions at these two sites.
Shiau-Chiou is a possible disposal site for low-level radioactive waste in Kinmen County. The groundwater quality data measured by the Industrial Technology Research Institute (Table 11) were used in this study; the data were entered into the four models in order to calculate the activity coefficients of the solution ions using the Davies and modified Debye-Hückel equations. The ionic strength calculated for the groundwater samples obtained from Shiau-Chiou ranged from 0.3 to 0.5 mol/kgw (Figure 17), indicating that Activity coefficient Using modified Debye-Hückel eq. for different ionic strength Figure 16. Dependence of the activity coefficient of Na + on the temperature calculated by the modified Debye-Hückel equation through EQ3/6. Table 11. Groundwater quality of Shiau-Chiou, Kinmen County, and Bai-Shi Lake, Hsinchu County.   6. Discussion

Internal Consistency of TDBs for Supporting the Safety Assessments of Radioactive Waste Repositories
There is no internal database for PHREEQC or any other geochemical model. The different geochemical models are shipped with different thermodynamic databases. Each database is derived from different data sets or historical collections of data and experiments that possibly focus on systems with different conditions (e.g., the ionic strength and the pH of a solution, the pressure state and the temperature range), and they include dif- The data of the water quality obtained at Bai-Shi Lake (Table 11) were compared with those obtained at Shiau-Chiou. The results showed that the water quality at Bai-Shi Lake complied with potable water standards, whereas the groundwater at Shiau-Chiou had higher concentrations of ions, with I ranging between 0.3 and 0.5 mol/kgw. The ionic strength levels calculated at Bai-Shi Lake by the four models ( Figure 19) were <10 −4 m, and the activity coefficient was approximately 1, suggesting that the concentrations did not need correction through the activity coefficient models.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 26 of 36 Figure 18. Activity coefficients of the major species for the four models with different activity coefficient equations at the Shiau-Chiou site. Figure 19. Calculated ionic strengths for the four models with different activity coefficient equations at the Bai-Shi Lake site.

Internal Consistency of TDBs for Supporting the Safety Assessments of Radioactive Waste Repositories
There is no internal database for PHREEQC or any other geochemical model. The different geochemical models are shipped with different thermodynamic databases. Each

Internal Consistency of TDBs for Supporting the Safety Assessments of Radioactive Waste Repositories
There is no internal database for PHREEQC or any other geochemical model. The different geochemical models are shipped with different thermodynamic databases. Each database is derived from different data sets or historical collections of data and experiments that possibly focus on systems with different conditions (e.g., the ionic strength and the pH of a solution, the pressure state and the temperature range), and they include different elements, minerals and ionic species. The chemical thermodynamics databases collect the values of the selected formations (Gibbs energy, enthalpy, entropy and constant pressure heat capacity), reaction data (equilibrium constant, Gibbs energy, enthalpy and reaction entropy), their uncertainty intervals, and the species and reactions in the safety assessments of radioactive waste repositories [55][56][57]. In addition, the modification of the activity coefficient by ionic strength must follow the activity coefficient equation, such as the Pitzer, Davies, modified Debye-Hückel equations, and the specific ion interaction theory (SIT) model. Thermodynamic data are collected from different data sets or experimental data. The internal consistency of a thermodynamic database must account for at least six different requirements: (1) all of the data in the database must conform to the basic thermodynamic relations and their results; (2) all of the data must be derived from a set of reference values (such as the reference temperature and pressure, the thermodynamic properties of elements and other basic compounds, etc.) and constants (such as the gas constant and molecular weight, etc.); (3) the standard state should be properly selected and suitable for all similar substances; (4) the appropriate mathematical model must be selected to fit all of the temperature and pressure data; (5) the appropriate aqueous chemical model must be selected to fit all of the aqueous solution data; and (6) all of the relevant experimental data and the data types must expressed and considered, and the conflicts and inconsistencies between the experimental measurements must be resolved [30,31].
The thermodynamic databases supporting the safety assessment of repositories may be internally inconsistent because they usually come from two or more source databases. Generally, the consistency of their internal databases can be verified through the following procedures: [31] • Check for internal consistency 1. Formation data The standard Gibbs energy of formation (∆ f G 0 X ), enthalpy of formation (∆ f H 0 X ) and entropy (∆ f S 0 X ) for the chemical compound X are in the standard condition (T = 298.15 K). The internal consistency of the formation data for the chemical compound X is checked with the Gibbs-Helmholtz equation: ∆ f S 0 X . is the absolute entropy of the chemical compound X in the standard condition, α i is the stoichiometric coefficient of the component element i in the chemical compound, S 0 i is the absolute entropy of the element i divided by the stoichiometric coefficient of at standard condition, z X is the charge of the chemical compound X in the standard condition, and S 0 H2(g) /2 is the electron entropy by convention.
For chemical compound X, the standard Gibbs energy (∆ f G 0 X ), enthalpy (∆ f H 0 X ), and entropy (∆ f S 0 X ) are calculated by Equation (16). All of these values must conform to Equation (15) in the standard condition.

Reaction data
The relationship between the standard Gibbs energy and the equilibrium constant must conform to Equation (5) at T = 298.15 K.

•
Check cross-consistency 1. The consistency between ∆ r G 0 X and the ∆ f G 0 X of the subject species can be expressed as: where the first and second terms are the total accumulation of all of the product species (P) and all of the reactants (r) for standard Gibbs energy, and υ P and υ r are the stoichiometric coefficients of the products (P) and reactants (r) involved in the reaction, respectively. The consistency in the relationship among ∆ f G 0 X , the product term ∆ f G 0 X,P and the reactant term ∆ f G 0 X,r can be determined by Equation (17). 2. The consistency between ∆ r S 0 X and the S 0 X,P of the subject species can be shown as: where the first and second terms are the total accumulation of all of the product species (P) and all of the reactants (r) for entropy, and υ P and υ r are the stoichiometric coefficients of the products (P) and reactants (r) involved in the reaction, respectively. The consistency in the relationship among ∆ r S 0 X , the product term S 0 X,P and the reactant term S 0 X,r can be determined by Equation (18).
Therefore, if the problem of accuracy and internal consistency in the thermodynamic data has been solved, the absolute accurate thermodynamic data of each solution species and mineral and the accurate geochemical model of the activity coefficient have been obtained. However, due to the complexity of real natural systems, predictions based on consistent thermodynamic data and accurate geochemical models may not be able to fully approach the repository system. Such a problem lies in the fact that geochemical models mostly rely on the assumption of local equilibrium. It is necessary to add the nonequilibrium kinetic reaction model and the data on mineral formation or the adsorption reaction, and even to couple the heat, pressure and stress fields into the model. Therefore, the development of geochemical models and the establishment of complete, consistent thermodynamic databases and kinetic reaction data still need continuous development in order to effectively carry out complete safety assessments of radioactive waste disposal site.

Using SCM to Estimate Kd for Performance Assessment
In the past, the traditional method for the determination of the retardation factor (R f ) of a nuclide released from a repository involved sampling the groundwater and geological materials of the disposal site. The samples were then used to carry out adsorption experiments and estimate the nuclide Kd. After comparing the experimentally derived Kd value with data from the scientific literature, a set of recommended Kd values were generated for the safety assessment of the disposal site. The Kd value is the ratio of the mass of the nuclide adsorbed or precipitated onto the minerals per unit of dry soil mass, S, to the concentration of the nuclide in the liquid, C. Equilibrium conditions and complete reversibility are assumed to exist between the soil minerals and the nuclides in aqueous solution. Kd can be expressed as Kd = S/C. R f is defined as R f = 1 + (ρ b /n e )Kd = v p /v c , where ρ b is the bulk density of the porous medium, n e is the effective porosity of the medium at saturation, v p is the velocity of the water traveling through the porous medium, and v c is the velocity of the contaminant traveling through the porous medium.
Many studies has pointed out that Kd values can be measured from many adsorption processes; however, experiments have shown that Kd values are not constant, and that they depend on the soil properties. As the geological material of the disposal system is a dynamic system, the Kd values will not remain constant over the long term. The soil properties affecting the Kd value include the soil texture (sand, loam, clay or organic soil), the soil organic matter content, the soil pH, the soil solution ratio, the solution or pore water concentration, and the presence of competing cations and adsorption sites for surface complexation. Because of the soil characteristics, the distribution coefficients of radionuclides in various soils may differ by several orders of magnitude. Table 12 lists the Kd values of nuclides in four representative soil types [58,59] and compares these values with the Kd values estimated in the case study (Table 6). Bulk MX-80 contains 65-75% montmorillonite, 10-14% quartz, 5-9% feldspar, 2-4% mica and chlorite, 3-5% carbonates and chlorite, and 1-3% heavy minerals. Therefore, montmorillonite and quartz are the main components and the dominant reactants in bentonite [60]. It can be seen from Table 12 that the variation of the Kd values is very large in sand, loam, clay and organic soils. Bentonite is mainly composed of montmorillonite clay and quartz sand; thus, the case study estimation in this study is between the Kd values of clay and sand proposed in the literature. Because the montmorillonite content of bentonite is 65-75%, the Kd value is estimated to be 28.6-33.0 cm 3 /g. Thus, the Kd value can be deduced by applying the chemical thermodynamic data of SCM and the chemical properties of minerals. Coupled surface complexation reactions in reactive chemical transport models can be used to simulate the spatiotemporal characteristics of Kd distribution caused by changing chemical conditions. In this study, the simulation of SCM was used to analyze the adsorption data measured in the laboratory, and a quantitative relationship between the measured Kd and the nuclide adsorption degree and environmental parameters (such as the pH, component concentration and density of the adsorption sites) was established. In this way, the simulation of SCM can provide scientific support for the measured Kd, and can increase the confidence of the Kd estimation. The potential uses of SCM for performance assessment and safety analysis in radioactive waste disposal include the following: 1. The application of supporting specific Kd estimation in the safety assessment calculation framework; the use of SCM simulation to directly estimate Kd can provide approximate values of parameter ranges. The interpretation and demonstration of the measured data indirectly supports the rationality of the experimentally derived Kd. As for the variation and uncertainty of the geochemical conditions, the sensitivity and uncertainty of the Kd values are analyzed, and the range of Kd is calculated.
2. Guiding the development of the experimental analysis method for obtaining Kd; using the chemical characteristics of surface complexation reactions, we can perform a sensitivity analysis of the environmental conditions for Kd. The experimental design was optimized by the Kd screening calculation in order to establish the most critical localization conditions of the Kd experimental measurement.
SCM simulation was applied in order to obtain Kd due to the importance of Kd in safety analysis and the influence of external environmental conditions on the Kd values. Due to the uncertainty and possible variability of the geochemical conditions associated with radioactive disposal sites, as well as the complexity of natural and engineering barriers, it may be impossible to obtain the Kd of nuclides by experimental methods for each region and all geochemical environmental conditions. Therefore, according to the experimental conditions of specific sites, more chemical data may be used to approximate, simplify, or generalize the data analysis to approach the site environmental state in order to establish confidence in the assessment of the radionuclide adsorption with different Kd values under various conditions through safety analysis cases in a larger spatial and temporal range.

Deterministic and Probabilistic Method for the Safety Assessment of Radwaste Disposal
The techniques for the evaluation of the evolution of nuclide geochemistry include two types of method: deterministic methods and probabilistic methods. The input parameters for deterministic and probabilistic methods should reflect the concentration of the adsorption or total dissolution of the nuclide of interest and the concentration of the released nuclides (e.g., high or low Kd values, the adsorption capacity, and the released concentration value). Because deterministic methods use point estimation to represent input parameters (e.g., chemical component concentration, equilibrium constant, temperature, pH and pe), using average values to represent these parameters may lead to the underestimation of the calculated nuclide concentrations, while the use of upper bound estimates would result in overly conservative estimates of the simulation scenarios. Meanwhile, probabilistic methods use the estimation of a probability density function (PDF) to estimate the scenario development profile, and can generate parameter distributions and uncertainties (e.g., the distribution of Kd and chemical equilibrium constants). Both the aleatory uncertainty and epistemic uncertainty associated with the long-term safety assessment of radioactive waste disposal often exist in engineering practice. However, the difficulty in dealing with the probability distributions of the parameters included in the long-term safety assessment of radioactive waste disposal can mean that there is no reasonable probabilistic safety assessment for radioactive waste disposal. A general method for safety assessment can be achieved through the statistical processing of the measured parameters. The aleatory uncertainty is quantified as a probability distribution. Additionally, based on expert judgment, the epistemic uncertainty is quantified as a probability distribution [61][62][63][64].
Crawford [62] proposed the modified transfer factor to establish the Kd value probability distribution function for the evaluation and application of Kd, including the following: (1) a surface area normalization transfer factor (f A ), (2) a mechanical damage transfer factor (f m ), (3) a cation exchange capacity transfer factor (f CEC ), and (4) a groundwater chemical transfer factor (f chem ). Because the first three transfer factors-f A , f m and f CECare defined as the ratios of variables that have either normal or log-normally distributed uncertainty ranges, the transfer factors themselves also have log-normally distributed uncertainty ranges. The term f chem is a groundwater chemistry transfer factor that relates the theoretical value of the distribution coefficient under the application conditions, K d(app), to the theoretical distribution coefficient value for a reference groundwater, K d (ref) . The f chem = K d(app) /K d (ref) . Because it contains uncertainty and spatially variable (usually nonrandom) variables, there is no reason to make it a lognormal distribution. However, f chem can be estimated from the functional relationship between the ionic strength and f chem , and then the convolution is performed between the probability density function of Kd and the probability density function of f chem . Thus, the uncertainty range of the conditional Kd value corresponding to the specific hydrochemical characteristics of groundwater is obtained [65].
With the development of safety assessment techniques for radioactive waste repositories, the demand for support for the modeling of adsorption data is increasing. In the future, it is expected that the experimental determination of radionuclide Kd will still be the main means for the estimation of the adsorption characteristics of the whole repository system. By using SCM chemical thermodynamic modeling to guide and establish the experimental determination of Kd, more parameter data and higher-confidence Kd values can be obtained. In the future, it is feasible to apply surface complexation reaction models to the study of the adsorption characteristics of radionuclides and the safety assessment of radioactive disposal sites.
The geochemical simulation applied in this study for the safety assessment of radioactive waste disposal is a deterministic model. The parameters input into the model are considered not to contain probabilistic components, and each component and input parameter is accurately determined by observation or experiment. However, due to the complexity of the disposal system, the possible inadequacy of in situ data collection, and the uncertainty introduced by the temporal variations in the system and its relevant parameters, a geochemical model is a simplified representation of the natural phenomena. In many cases, probabilistic models are used to simulate deterministic systems, and deterministic models are used as subunits in a larger stochastic framework, including small-scale phenomena that cannot be observed or simulated accurately.
Although the traditional deterministic method for the safety assessment of radioactive waste disposal can ensure sufficient safety, the probabilistic method can be used to quantify the performance assessment more clearly. In order to improve the quantitative evaluation of uncertainty in the performance assessment of radioactive waste disposal, researchers should consider the package in terms of deterministic and probabilistic performance assessment, including what could go wrong, how likely this occurrence is, and what the consequences will be. The combination of deterministic and probabilistic methods can help us to evaluate the performance and risk possibility of a radioactive waste disposal site, and can also be applied to decision-making [66].

Application of Geochemical Modeling for the Performance Assessment of Radioactive Wase Disposal
Knowledge on the species distribution of radionuclides is required in order to understand the state equilibrium and geochemical reactions at a disposal site. The radionuclide compositions and groundwater quality data can form the basis for reactive chemical transport simulation for the determination of the transport mechanisms and migration phenomena of radionuclides. In order to better represent and predict the chemical reaction and behavior of radionuclides in radioactive waste repositories, laboratory tests, field experimental studies and geochemical simulation studies should be drawn upon. Despite the significant improvements made in geochemical models, there is still potential for further development.
Before the execution of a geochemical model simulation, the following items should be established: computer code, a database of thermodynamic data, the kinetics of chemical species, and measured data of the chemical and physical properties in the study area. Accordingly, this study proposes that the following questions should be considered: (1) Are the chemical analysis data sufficiently accurate to support simulation studies? (2) Does the TDB used in the program include data on the key solution species and reactions of surface complexation models, which occur between chemical species and minerals? (3) Are the equilibrium constant values of crucial reactions in the TDB sufficiently accurate to support a clear and unambiguous simulation conclusion? (4) Is the ionic strength level sufficiently low to ensure the appropriate applicability of the activity coefficient model in the geochemical program?
On the basis of the above questions, the following recommendations are proposed for future research.
(1) Before establishing a detailed geochemical model, researchers must determine which model to develop; because most computer models are based on the assumption of chemical equilibrium, it is important that the assumption of chemical reaction equilibrium (in disposed radioactive waste) should be evaluated and confirmed using kinetic simulation models. The time required for a reaction to reach equilibrium can be calculated using a kinetic model; subsequently, it can be decided whether the assumption can be justified.
The problem must be thoroughly understood before the simulation, followed by suitable simplification (if possible) to create an appropriate conceptual and geochemical model. The appropriate simplification of the conceptual simulation not only reduces the problem's complexity considerably but also enables the achievement of the desired results. Therefore, before the establishment of a geochemical model, the following points should be considered: (a) the species and components to be included in the model; (b) the types of geochemical reactions to be included in the model; (c) whether chemical reaction kinetics are included in the model or a local equilibrium hypothesis is used; and (d) whether any other processes need to be included in the simulation.
(2) Activity coefficient calculation: the selection of a suitable chemical TDB is crucial in the generation of an appropriate species activity coefficient. The ion-specific parameters provided in a TDB are highly crucial for the applicability of the modified Debye-Hückel equation. Because thermodynamic data affect the species distribution pattern generated in the model simulation, the simulation results will reflect the accuracy of the activity coefficient. For example, in a dilute solution (I < 0.1 mol/kgw), if the TDB is not accurate, the calculated ionic strength will deviate. This may lead to erroneous activity coefficients of high-charge species. Thus, the careful selection of TDBs is necessary.
Pitzer developed a theoretical model of an electrolyte solution based on the Debye-Hückel equation, which is highly effective in calculating the activity coefficients of highconcentration mixed electrolyte solutions. The model does not contain parameters to modify the standard state or consider independent reactions among chemical species. Its application is currently limited to low-temperature and low-pressure conditions. The Pitzer equation has been used in several geochemical simulation software packages. It is often used for the analysis of environmental problems involving high-concentration solutions.
(3) The thermodynamic equilibrium constants of a reaction must be further supplemented and improved. Currently, the equilibrium constants of various minerals and soils for surface complexation models of radionuclides are scant. Hence, experts should compile and evaluate data in order to determine the recommended values of such constants. The composition of soil is complex and varies according to the place of origin, which imparts strong regional characteristics to the equilibrium constants of surface complexation models. Therefore, it is likely that researchers will need to customize their databases to their local conditions.
(4) Surface complexation should be included in chemical transport modeling in order to evaluate the safety of nuclear waste disposal. Empirical methods, such as the isothermal linear model (Kd method) [67][68][69][70], were previously used to calculate the groundwater and radionuclide transport in the performance assessment of radioactive waste disposal. Although the Kd method can describe the heterogeneous interface behavior of radioactive species at specific pH and Eh, as well as the composition of the ions in the surrounding water, it is not suitable for situations involving changing geochemical conditions. On the other hand, a surface complexation model can account for changes in hydrochemistry and mineral ion composition with time, and thus better represents the effects on migration of radionuclides under varying geochemical conditions [71]. Surface complexation should be included in chemical transport models in order to better support the reactive chemical transport of radionuclides between water and minerals. Thus, it provides a more powerful simulation basis for heterogeneous interfaces in geochemical reactions.

Conclusions
Maintaining the low solubility of radionuclides in the geochemical environment of groundwater systems is an effective way to retard the release of radionuclides. In order to better understand and improve the performance evaluation of nuclear waste disposal facilities, quality thermodynamic data on water-mineral and nuclide geochemical reactions in geochemical environments are essential to ensure a strong scientific basis for the evaluation of the safety performance of nuclear waste management. Combining deterministic and probabilistic methods can help us to evaluate the safety and risks associated with a repository's performance. This information can then be applied to the decision-making process for radioactive waste disposal.
The site selection, construction, monitoring and management for radioactive waste disposal are complex processes that require a multidisciplinary approach involving geoscientists, engineers, environmental scientists, socio-economists and cultural impact assessment experts. The stakeholder response in complex decision-making processes for the performance assessment of radioactive waste disposal facilities and nuclear waste management is continuously being developed. The results of this study can be used as technical discussion issues of stakeholder participation in the decision-making process of radioactive waste disposal.
This study reviewed the research and application of previous geochemical models, including the development of geochemical models, the source of TDBs and their application in geochemical code, the adequacy of TDBs for surface complexation models and case studies, and the selection and application of activity coefficient models in geochemical models. In addition, case studies were conducted and differences in the activity coefficients between geochemical models were examined. The PHREEQC, MINEQL+, MINTEQA2, and EQ3/6 models were used in this study. The study results demonstrate that when the solution ionic strength was <0.5 mol/kgw, the differences between the activity coefficients calculated using the Davies and modified Debye-Hückel equations were <5%. The difference between the Pitzer and Davies equations, or those between the Pitzer and modified Debye-Hückel equations in terms of the calculated activity coefficients were <8%. The effect of the temperature on the activity coefficient is fairly small over the range of temperatures considered (i.e., 20°C to 90°C), and only slightly influenced the model results of the Davies and modified Debye-Hückel equations. Two field cases were also examined. The ionic strength was determined to be <0.5 mol/kgw for the Shiau-Chiou groundwater, such that the Davies or modified Debye-Hückel equation can provide accurate activity coefficients for such water. The ionic strength of the water from Bai-Shi Lake was calculated to be <10 −4 mol/kgw, suggesting that the solution activity coefficient was close to 1 and did not require correction.
In this study, we reviewed the techniques and methods used in geochemical model-ing in relation to nuclear waste disposal, and we established a process for the estimation of Kd using SCM. Both deterministic and probabilistic models can be used to analyze the safety evaluation of repositories. In the future, the probability distribution and uncertainty of the parameters of Kd and equilibrium constants can be used in geochemical and reactive transport models to simulate the long-term safety of nuclear waste disposal sites. The study results can provide a strong scientific basis for the assessment of the safety of nuclear waste repositories, and support the development of environmental management or remediation schemes for sites with near-surface contamination.  Institutional Review Board Statement: The study did not involve humans or animals.

Informed Consent Statement:
The study did not involve humans.