Mathematical Modeling of a Cs ( I ) – Sr ( II ) – Bentonite – Magnetite Sorption System , Simulating the Processes Taking Place in a Deep Geological Repository

The derivation of mathematical models of systems consisting of Cs(I) or Sr(II) and of bentonite (B), magnetite (M) or their mixtures (B+M) are described. The paper deals especially with modeling of the protonation and sorption processes occurring on the functional groups of the solid phase, namely on so called edge sites and layer sites. The two types of sites have different properties and, as a result, three types of Surface Complexation Models (SCM) are used for edge sites, viz. two electrostatic SCMs: the Constant Capacitance Model (CCM) and the Diffusion Double Layer Model (DLM), and one without electrostatic correction: the Chemical Model (CEM). The processes taking place on the layer sites are described by means of the Ion Exchange Model (IExM). In the course of modeling, the speciation of the given metal in the liquid (aqueous) phase has to be taken into account. In principle, the model of protonation or sorption processes is based on the reactions occurring in the aqueous phase and on the surface of the solid phase, and comprises not only the equations of the equilibrium constants of the individual reactions, but also the mass and charge balance equations. The algorithm of the numerical solution is compatible with FAMULUS 3.5 (a Czech software product quite extensively used at Czech universities in the last decade, the bookcase codes of which are utilized).


Introduction
In our deep geological repository, we plan to use canisters made of stainless and carbon steel, with a compacted bentonite barrier.The use of bentonite as a backfill barrier in repositories for nuclear waste is based mainly on its low permeability, swelling properties and capability to significantly retard the migration of most radionuclides (RN) [1,2,3].As a result of corrosion processes, the main corrosion product of steel canisters -magnetite -will also form an important part of the barrier.To be able to predict the fate of RN in repositories, the retardation processes, i.e. mainly the sorption processes, have to be studied, and mathematical sorption models of the 'RNbentonite -magnetite' system have to be developed.
The bentonite (clay) surface contains at least two types of surface groups.The first type are permanently charged functional groups created by ionic substitution within the crystal structure.Isomorphic substitution of, e.g.Al 3+ for Si 4+ within the tetrahedral layer creates a permanent negative charge on the mineral surface, which is compensated externally by cations.These inter-structural-charge surface sites are denoted as layer sites.On the edges of the surface structure, there are (ºSOH) sites with a pH-dependent charge.This is due to the "adsorption" of H + ions (then so called protonation proceeds: ºSOH ® ºSOH 2 + ) or "desorption" of H + ions (then deprotonation proceeds: ºSOH ® ºSO -) depending on the pH of the solution.These variable-charge surface sites are designated as edge sites.These two surface site types are responsible for two uptake processes.The first, on layer sites, tends to be dominant at low pH and/or high sorbate concentrations.The mechanism of this process is cation exchange.
The second process occurs on edge sites, depending on pH, the mechanism of which is surface complexation.As for magnetite, both types of surface sites were also found [4].The aim of this work was to derive the mathematical models of the protonation and sorption processes occurring in the 'RN -bentonite -magnetite' system using the FAMULUS software product, including the Newton-Raphson nonlinear regression method [5], and to prepare codes for evaluating the experimental data.

Types of models and basic modeling approaches
In principle, four types of sorption models can be used.The simplest is the K D -model characterizing the linear equilibrium isotherm.The second, e.g. the Langmuir equation, describes the non-linear sorption isotherm.The third is based on the application of a classical chemical-equilibrium equation or equations, e.g., the Ion Exchange Model (IExM).The fourth type is represented by Surface Complexation Models (SCMs) [6,7,8].As was mentioned above, IExM and SCMs describe the processes occurring on layer sites and edge sites, respectively, and as a result, they can be designated as the most important models.
At least three types of surface complexation models (SCMs), namely two electrostatic models, the Constant Capacitance Model (CCM) and the Diffuse Double Layer Model (DLM) and one chemical equilibrium, non-electrostatic model (CEM), were employed to simulate the amphoteric property of the solid phase types of bentonite, goethite and montmorillonite, and to describe the sorption of different metal ions and their complex compounds from an aqueous solution as a function of pH, ionic strength, solution concentration, etc. Well developed SCMs enable the behaviour of systems, e.g.RN-bentonite-magnetite, to be quantitatively described and the sorption processes to be simulated.
SCMs and IExMs are based on the following suppositions: l No mutual interactions exist between particles adsorbed on the surface of the solid phase.
l The protonation and deprotonation processes depend on pH.
l On the surface of the solid phase (e.g.bentonite), as was mentioned above, there are two types of functional groups, i.e. so called edge sites and layer sites.
l The functionality of adsorption sites, i.e. edge sites or layer sites, is identical.
l The concentration of the ith-component in the aqueous phase near the surface (in the aqueous layer adhering to the surface) is given by the Boltzman equation ( 1): where (C i ) S is the concentration near the surface, C i is the so-called bulk concentration of the ith-component, z is the charge of the ith-component, y is the electrostatic potential, F is the Faraday constant, R is the gas constant and T is the absolute temperature.
l Between the surface charge, s, and the electrostatic potential, y, the following dependencies hold (2,3,4): where G is the so-called Helmholtz capacitance.
In order to apply SCMs for describing the sorption processes occurring in the system type of e.g.RN-bentonite--magnetite, we need the surface area of each solid phase component, the total surface site concentrations (edge and layer sites) and the protonation constant (for edge sites) and ion exchange constant (for layer sites).These can be obtained by evaluating the acid-base titration data (in our case of bentonite and magnetite and their mixtures) and by measuring the surface area.In the course of evaluating the of sorption experimental data, specifically the sorption dependencies on pH, the values of the equilibrium constants of the surface complexation reactions are found.
As for the modeling approaches used in the case of RN--bentonite-magnetite type system, two different procedures can be tested: l The Generalized Composite Approach (GC), where the given mixture of solid phase components, e.g.bentonite + magnetite, is considered as a compact sorbent characterized by a single set of titration and sorption parameters, which are sought by direct fitting of the experimental data.
l The Component Additivity Approach (CA), composed of a weighted combination of models describing the protona-tion + ion exchange (on layer sites) and the sorption on individual solid phase components, e.g. on bentonite and magnetite.The individual parameters have to be obtained by fitting the appropriate experimental data that are valid for the solid component.
It is evident that the GC Approach demands less experimental time than the CA Approach because we need only one titration curve and one sorption pH-dependency for the given mixture of solid components.On the other hand, if the protonation + ion exchange and sorption quantities characterizing the individual solid phase (mineralogical) components are determined, the CA Approach enables the sorption behaviour of selected component mixtures to be simulated.It is also evident that the codes corresponding to both the GC and the CA approaches are in principle different: whereas the GC Approach code has to be based on a non-linear regression procedure, the CA Approach code is relatively simple, even if internal iterative loops also have to be used.

Description of the system
where SO -, SOH and SOH 2 + are symbols for edge sites, and X -is the symbol for layer sites.
In the course of titration, protonation reactions ( 5) and ( 6) on the edge sites, and the ion exchange reaction (7) on the layer sites are under way in the system.Titration occurs under an inert atmosphere (e.g.N 2 ), and as a result, the presence of atmospheric CO 2 need not be taken into account.The determination of the titration curve starts at approx.pH 7, namely the titration of an aqueous suspension of the solid phase, having a given ionic strength, with both HCl and NaOH solutions such that the total titration curve consists of two parts -two sets of experimental points obtained with both HCl and NaOH.

Derivation of CEM, CCM and DLM titration models
It is supposed that the electroneutrality represented by Eq. ( 8) exists between positive and negative charged surface groups and negative and positive charged species in solution in the course of titration.
By rearrangement of Eq. ( 8) we can obtain Eq. ( 9).It describes the course of titration relating the surface charge, Q [mol/kg], and the experimentally measurable parameters: where Eq. ( 9) consists of two parts: the right-hand side can be designated as the experimentally determined values of the surface charge, Q exp , and the left-hand side expresses the sum of the two values of the surface charges, and the charge of the layer sites, , np, describes the experimental titration curve, consisting of np experimental points.The goal of modeling the titration curve is to construct the function ) applicable for fitting of experimental data and for determining of the values of the protonation and ion exchange constants, and the total concentrations of the edge sites and layer sites in the given solid phase.As for the numerical method, the acid-base titration data are fitted by the Newton-Raphson multidimensional non-linear regression method and the quantity WSOS/DF (see Eqs. ( 83) and ( 84)) is used as the criterion for the goodness-of-fit.

Derivation of Q cal = f([H + ]) using CEM (Chemical Equilibrium Model)
As was mentioned above, it is assumed that in the case of CEM y = 0 and as a result, the concentrations of the ith--components existing near the surface, (C i ) S , equals the bulk concentrations, C i (see Eq. ( 1)).
The protonation constants of the edge sites, KS1 and KS2, are given by Eqs. ( 12) and ( 13), respectively.Now, using the balance equation ( 14) for the total concentration of the edge sites, SSOH, together with Eqs. ( 12) and ( 13), we derive Eq. (15) corresponding to the function The surface charge can also be expressed in coulomb per m 2 , s ES [C×m -2 ]: As regards layer sites, the ion exchange reaction takes place on these sites, and the equilibrium constant, KS5, is given by Eq. ( 17).Furthermore, the dissociation of XH and XNa according to equations ( 18) and ( 19), respectively, needs to be taken into consideration.On the basis of the literature data, it is supposed that for the corresponding values of the dissociation constants, it holds: KS5b KS5a and KS5a £ 10 -2 .It follows from this that the dissociation of XH can be neglected and that the concentration of XNa is very low and therefore XNa is practically dissociated.
If the function ) is to be derived, two balance equations are needed, viz., the first balances the layer sites, cf.Eq. ( 20) for total concentration of layer sites SX, and the second balances the sodium ions, cf.Eq. ( 21) from which Eq. ( 22) can be obtained. [ ) In our case, the starting concentration [Na 0 ] equals the starting value of the ionic strength achieved by adding NaNO 3 , or NaClO 4 and as for [XNa 0 ], if the starting pH of titration is approx.7, it is assumed to be approximately equal to the value of SX (cf.Eq. ( 20)).
Also in this case, the surface charge can be expressed in Coulomb per m 2 , s LS : It is evident that the function ) is given by the combination of Eqs. ( 23) and ( 22).
Altogether, the CEM function of the titration curve, namely ), which can be used for evaluating the experimental data, consists of equations (15), ( 22) and (23).Formally, the total surface charge in Coulomb per m 2 , s S , can be given by Eq. (25): In principle, the construction of this model function is congruent with the construction of the CEM model.However, the value of electrostatic potential y does not equal zero and Czech Technical University in Prague Acta Polytechnica Vol. 45 No. 5/2005 the Boltzman equation, Eq. ( 1), has to be used for calculating the component concentrations existing near the surface, (C i ) S .These are then inserted into the model equations ( 15), ( 22) and (23) instead of the bulk concentrations.First, quantity y must be calculated, namely by means of equations ( 2) and (25).After the rearrangement procedure, equation ( 26) is obtained, for the solution of which a suitable interpolation method has to be used.
[ ] [ ] where the symbols with subscript "s" mean the concentrations near the surface.The proper CCM model function consists of four equations, namely (26), ( 27), (28a) and (28b).The values of the quantities searched, viz.KS1, KS2, KS3, SSOH, SX and G, are determined in the course of simultaneous solution of the above mentioned equations using a non-linear regression procedure.

Description of the system
The derivation of SCMs (CEM, CCM, DLM) for sorption of strontium on bentonite or magnetite will be demonstrated, as an example.The sorption experiments are carried out in a mixed batch reactor under given conditions, i.e., under given values of the starting volume, V 0 , and composition of the aqueous phase (as a starting solution, synthetic granitic water is used) and mass of solid phase, m.The value of the given ionic strength, I, is adjusted using NaNO 3 .The reaction time (time of contact of the phases, e.g.approx.30 days) must be sufficient for equilibrium to be reached.Because the system is open to the atmosphere, the influence of atmospheric partial pressure of CO 2 is taken into consideration, especially if pH is greater than approx.8.2.In the course of the sorption experi-ment, the influence of pH is observed, namely in such a way that each experimental point has a given pH, the value of which is adjusted by means of 0.1 M HCl or 0.1 M NaOH.
The sorption phenomenon is observed to proceed with the formation of a surface complex, or complexes, including ion exchange on the layer sites, as described below, between the surface groups and various species of strontium (Sr 2+ , SrCO 3 0 , SrNO 3 + ) present in the experimental solution.These species compete with each other to form a surface complex with the solid phase, and the values of the corresponding equilibrium constants quantify this competition.The input data include among others the protonation and ion exchange constants and the total concentrations of the edge and layer sites, or the Helmholtz capacitance, determined in the course of evaluating the titration curve (KS1, KS2, KS5, SSOH, SX, or G).

Derivation of CEM, CCM and DLM sorption models
The experimental sorption data are in the form: Kd exp = f(pH exp ) or %Sorption exp = f(pH exp ), where Kd exp is the distribution coefficient of sorption of Sr(II) and %Sorption exp expresses the sorption of Sr(II) in percentage units.As a result, the analogous model function, namely Kd cal = f(pH) or %Sorption cal = f(pH), has to be derived using the following procedure.
Firstly, let the reactions occurring on the solid phase (e.g.bentonite) and in the aqueous solution be formulated (the symbols for the corresponding equilibrium constants are given in parenthesis): Secondly, the equilibrium constants of the reactions mentioned above have to be established: Constants holding for edge-sites (using the Boltzman equation (1)): Constants holding for layer sites: Constants of reactions occurring in an aqueous solution (for Where SR is the solubility product of SrCO 3 .solid and pCO 2 is the partial pressure of CO 2 (in our case, it deals with atmospheric partial pressure).
The Davies equation ( 64) for the i-th ionic species and the Hegelson equation ( 65) for the i-th neutral species are used for calculating the activity coefficients by means of which the values of the equilibrium constants are corrected for the given ionic strength, I.
where A = 0.509 (holds for aqueous solutions and ambient temperature), log f b I i = × (where b = 0.01-0.10). (65) Thirdly, the balance equations are formulated: Balance equation of the surface charge, s: It is evident that only charged surface species are taken into consideration.As was described above in connection with the derivation of the CCM and DLM models, it is necessary the quantity y needs to be determined, namely by solving Eqs. ( 67) and (68) originating from combination of Eq. (66) with Eq. ( 3) or (2), respectively.Depending on the type of surface complexation model, the obtained value of y is then inserted into Eqs.( 46 Balance equation of nitrates: Balance equation of sulphates: It is assumed that sulphates are present only in solution.
Balance equation of carbonates:

H CO SrCO
The solution of the set of equations mentioned above, namely for the given interval of pH, lies in combinating the balance equations with the equilibrium constant equations, for example, the concentration [SOSr + ] derived from Eq. ( 50) is inserted into Eqs.(66), ( 70) and (73), [XSrNO 3 ] from Eq. (55) into Eq.( 73) and (78), etc.The group of thus modified balance equations creates the regression function, the solution algorithm of which is depicted in Fig. 1.
It is evident that the algorithm consists of two loops: external and inner.In the external loop, the Newton-Raphson multidimensional non-linear regression procedure is used for fitting the experimentally determined data, and in the inner loop, the proper regression function is solved.In the last step, the functions Kd cal = f(pH), Eq. ( 79), and %Sorption = f(pH), Eq. ( 80 Sr SrNO Sr 0 100. (80) The sorption on the edge sites, Eq. ( 81), and layer sites Eq. ( 82), can also be calculated The goodness-of-fit is evaluated by the c 2 -test, which is based on calculating the quantity c 2 according to Eq. (83): where (SSx) i is the ith-square of the deviation of the experimental value from the calculated value, (s q ) i is the relative standard deviation of the i-th experimental point, N (º n p ) is the number of experimental points.The value of c 2 is used for calculating the criterion WSOS/DF (weighted sum of squares divided by degrees of freedom) [9]: where n i is the number of degrees of freedom, n p is the number of experimental points and n is the number of model parameters sought during the regression procedure.
If WSOS/DF £ 20, then there is a good agreement between the experimental and calculated data, while, in the other case the fitting is regarded as unsatisfactory.As for (s q ) i , on the basis of the corresponding experiments, the constant value for each experimental data point was assumed to be equal to 0.1 (in the case of the titration experiments) or 0.05 (in the case of the sorption experiments).
The codes were successfully verified using the two systems, namely Sr-bentonite-magnetite and Cs-bentonite-magnetite, and the results were published at the international conference MIGRATION'03 [10], then at two home conferences [11,12].The final paper has been prepared for MIGRATION'05 [13].
In principle, the same procedure and algorithm as was described above can be used for modeling uptake processes of the SCMs+IExM type that occur on hydrated metal oxides, clay materials, etc.
is the total volume and V 0 is the starting volume of the aqueous phase [dm 3 ]; m is the mass of the solid phase [kg]; [C a ] and [C b ] are the bulk concentrations of acid (i.e. of [Cl -]) and caustic soda (i.e. of [Na + ]), respectively, in solution [mol×dm -3 ]; c H and c OH are the concentrations, [mol×dm -3 ], of hydrochloric acid, HCl, and caustic soda, NaOH, respectively, used in titration; v H and v OH are the consumptions, [dm 3 ], of acid and caustic soda, respectively, in the course of titration.
where V is the volume of the aqueous phase, m is the mass of the solid phase, e.g.bentonite, and [Sr] 0 , [NO 3 -] 0 , [SO 4 2-] 0 and [CO 32-] 0 are the starting concentrations in the aqueous phase.If the partial pressure of CO 2 , pCO 2 , is taken into consideration then it holds (KpI is the constant of equation (63) corrected for the given ionic strength, I): by means of Eqs.(61) and (62), respectively.Balance equation of the total concentration of the layer sites, SX: