1.3.3　关键放射性核素在水岩体系中的种态及其分布

1.3.3　关键放射性核素在水岩体系中的种态及其分布

前已述及，放射性核素在不同的水岩体系中的存在种态及其分布是随体系的组成变化而变化的。因此，理解放射性核素在不同条件下的吸附、扩散、迁移等特性的前提是了解其种态存在形式及其变化。

核素的化学种态系指核素在特定水文地质条件下的具体存在形式，如和即为铀在某种条件下的2个种态。种态分析系指根据体系的组成、温度和压力，计算体系处于平衡状态时元素各种态的组成，即体系中各元素或核素的种态分布。自20世纪80 年代起，国外已经开发了多款化学种态分析软件，如PHREEQC、EQ3/6、MINTEQA2等。这些软件开发的先决条件是得到了在相应条件下的各种元素或核素种态的热力学参数，因此具有一定的地域性。由于在软件的编写过程中需要采用大量的热力学数据，因此，这些软件的编写工作首先是在美国等西方发达国家开始的。20世纪80年代，我国的一些科研人员在美国进修期间首次接触到这些软件，如EQ3/6，并在国内的一些工作中使用了这些软件。一般来讲，使用这些软件是相对比较简单的，只要在计算机上安装上相应的启动程序，并将程序要求的数据按一定的格式输入，便会得到结果。但是，对计算结果的解释往往存在很大的不确定性，有时甚至会出现一些预想不到的结果。与此同时，非软件编写国家的用户虽然可以使用这些程序，但往往不能得到源程序。因此，其计算结果只能是假定所计算的条件与程序中某些给定条件一致情况下的结果。毫无疑问，这多少是有点问题的。这也正是为什么国际上出现了不同款的种态分析软件。

我国是一个具有较宽广区域的国家之一，东西南北不同区域的水文地质条件有较大的差异，因此有必要开发适合我国特定条件的化学种态分析软件，以服务于我国核设施退役治理研发工作的需要，尤其是高放废物地质处置研发工作的需要。为此，北京大学核环境化学课题组在王祥云教授的鼎力帮助下，于2008年开始编写化学种态分析软件CHEMSPEC[12]。这是一项非常艰巨的工作。由于软件的编写需要具备一定的数学基础和软件编写能力，同时还需要对有关的化学反应有深刻的理解。王祥云教授是中国放射化学领域少有的既有坚实的化学基础，又有坚实的数学功底，且可进行程序编写工作的前辈之一。由于在程序的编写过程中，需要对实际体系进行计算，以便及时发现程序在运行过程中存在的问题并加以修正。课题组的陈涛等同学参与了这一工作。在王祥云教授的鼎力帮助下，北京大学核环境化学课题组编写了中国首款化学种态分析软件CHEMSPEC。这款软件经多名学生的使用，如陈涛、孙茂等利用该软件计算了Am 和Np在特定地下水中的种态分布和溶解度，其溶解度的计算结果与实验结果符合良好，证明了软件的可靠性[13-14]；朱建波等利用C++语言对软件界面进行了可视化，并增加了作图功能[15]；蒋美玲、蒋京呈、周万强等利用CHEMSPEC 计算了Am、Np、Pu、U在北山地下水中的种态分布，并利用新加入的表面配合模型模拟计算了关键核素在固液界面的吸附[1617]；兰图等利用CHEMSPEC计算了低浓缩铀靶辐照后溶液中铀的化学种态及主要裂变元素对铀种态分布的影响[18]。这在一定程度上对该软件的适应性进行了验证。

需要特别指出的是，这款软件所使用的数据主要来自国际上可公开获取的热力学数据库中的数据，与此同时，我们对一些可疑的数据进行了详细的分析和甄别。尽管如此，仍然有些数据需要进一步核实。

CHEMSPEC的基本原理是根据输入文件，基于热力学平衡调用计算程序并在数据库中搜索可能形成的所有种态，最终输出计算结果。在整个流程中，体系组成、平衡方程求解、反应模型和数据库是软件计算模型的关键组成。CHEMSPEC 的计算输入体系分为两大部分。第一部分是研究体系的环境条件，第二部分是研究体系物质的数量、存在形式、浓度以及是否参与氧化还原反应等信息。CHEMSPEC 将物质划分成种态和组分两类，化学元素可能存在的各种形式统称为种态，将可定量表示体系组成的最少种态称为组分。以Al（NO3）3-UO2（NO3）2-NaF-H2O 体系为例，可识别的输入组分为Al3+、F-、H+和H2O；以Am（NO3）3-Na2CO3溶液体系为例，体系的种态数目可多达19种，但考虑到各种态浓度同时满足多个平衡的关系，例如H2CO3和的浓度与电离平衡常数、H+和OH-浓度与水的离子积有关，最终完整描述该体系只需要6个组分。

正确读取给定的体系后，软件会进行种态求解。通过热力学数据计算体系平衡状态的两种方法分别是吉布斯自由能最小化法和质量平衡法。CHEMSPEC 程序选择利用质量平衡法求解化学反应方程。方程的具体构建过程如下：对任意组成为AaBbCc…的种态Sj，设其生成反应对应的热力学平衡常数Kj的表达式为式（1-18）。

式（1-16）～式（1-18）中a、b、c…为组分A、B、C 的化学计量因子。种态Sj的浓度引入活度系数γ 后的表达式为式（1-19），而每一个组分的质量平衡方程式为式（1-20）。其中，Ci为组分i的分析浓度，[Xi]为组分i的游离种态浓度，Xji为组分i在种态j中的化学计量因子，若体系有n 个组分，可以写出n 个质量平衡方程，联立求解可求出n 个游离组分的浓度，进而得到各种态的浓度。

lg K 与温度的函数关系用式（1-21）表示，A、B、C、D 和E为各种态的常数，程序中有25℃时的ΔH 值，其他温度下的lg K 由范德霍夫公式计算。

由一个种态的浓度计算其活度，需要知道其活度系数。活度系数的求算与体系的离子强度有关，根据不同的离子强度，程序会选择不同的公式进行计算。CHEMSPEC采用Debye-Hückel公式、Devies公式、B-dot公式和专属离子相互作用理论SIT公式和Pizer公式进行计算。国际上常用的种态分析软件EQ3/6与WATEQ 采取了不同的B-dot公式，Pizer公式适用于离子强度较大的情况，但相关参数极其有限。

对于中性物质的活度系数，CHEMSPEC 采用两种方法处理，一种是采用CO2的活度系数作为除H2O、H2S和Si（OH）4以外的中性溶液中其他各种态的活度系数，因为这4种物质在不同浓度的NaCl溶液中都有精确的实验值；另一种则是使用Setchénow 方程（1-22）：

计算H2O 的活度系数，采用适用于稀溶液的拉乌尔定律计算α 值。

式中，Waq为H2O 的质量，当Waq=1 kg时：

此外，规定难溶化合物（固态）、离子交换态、表面组分、静电组分及表面种态的活度均为1。

质量平衡方程的求解。求解质量平衡方程的本质是求解高阶非线性方程组，为达到快速、可靠收敛到有意义的解，CHEMSPEC 可在五种解非线性方程组的算法间自动跳转。首先调用牛顿-拉夫森迭代法，在多数情况下，都能快速收敛到合理解。如果不收敛，程序自动调用阻尼最小二乘法。在处理沉淀溶解平衡及氧化还原平衡问题时，难溶化合物的生成或价态变化往往使有些组分的游离浓度变得非常低，导致数值之间的数量级差异过大，不利于质量平衡方程组的求解，而阻尼最小二乘法可以避免迭代求解过程中的无效修正。如果仍不能收敛，程序将调用单纯型算法，这是一种无需求导的直接型最优化方法，虽然收敛得较慢，但对于初值的选取比较宽松。如果三种方法都不能收敛，程序调用广义逆法，以克服在计算一阶导数时出现的奇异问题。如果四种方法都不能收敛，程序将调用另一个直接搜索法——模式搜索法（或称步长加速法），轮番逐一进行一维搜索下降方向，若仍不收敛，迭代过程将被中断，此时程序显示迭代失败，询问用户是否要修改初值重新计算，或者就此终止。用户可尝试修改输入文件，或者变更被剔除种态的名单，以期达到收敛的目的。

根据数学求解，以各反应模型中的热力学平衡关系为桥梁，可以从实际体系中提取出计算方程并进行求解。目前，CHEMSPEC 能处理溶液体系中发生的酸碱反应、沉淀溶解反应、氧化还原反应、配位反应以及表面配位吸附反应。

（1）沉淀溶解平衡。具体地说，在处理沉淀反应中引入了饱和指数SI的概念，见式（1-25）。当SI＞0时，溶液处于过饱和状态；SI=0时，溶液恰好饱和；SI＜0时，溶液处于不饱和状态。因为事先无法知道何种物质将发生沉淀及沉淀的量，需要通过预报-校正或者迭代计算，使得最终的计算结果小于零或者等于零。CHEMSPEC 程序给出的最终结果是沉淀种类和具体的量，对于计算溶解度十分方便。

（2）氧化还原反应。在氧化还原反应中，体系的平衡关系与其氧化还原电位密切相关。对于任意氧化还原反应式（1-26），M 为氧化还原电对的主组分，S为从组分，m、n、h 和w 为化学计量因子，a 表示活度。设该反应的平衡常数为K，由热力学平衡可以得到其活度关系式（1-27）。

氧化还原电位可以由三种方式给出。一是直接给出体系的电位值或自由电子的活度，直接得到pe，如式（1-28）所示；二是由氧化还原优势电对来决定体系的氧化还原电位，反应满足浓度平衡关系式（1-29），进行转化后得到式（1-30）；三是由氧分压来确定，由反应方程可以得到pe的计算公式（1-31）。

在既有氧化还原反应，又有沉淀生成的情况下，程序采取“冻结”和“解冻”氧化态的策略。该策略首先将氧化态解冻，即调用子程序KERN 计算溶液的平衡组成时，容许所有可能的氧化还原反应发生，然后冻结氧化态，调用DEPOSI计算沉淀-溶解平衡，再解冻氧化态，调用KERN 计算溶液中包括氧化还原反应的平衡浓度，再次调用DEPOSI，如此反复，直至体系达到热力学平衡。这样就避免了同时计算氧化还原和沉淀反应所带来的复杂性，但是当体系中含有发生氧化还原反应产生的沉淀（RDXSS）时，由于RDXSS是伴随氧化还原反应而生成的，它既不能在“解冻”阶段，也不能在“冻结”阶段生成。程序为此采用了一种专门处理RDXSS的计算策略。

（3）离子交换模型。目前可用于解决比较简单的离子交换计算。假设初始吸附离子为+1价，可交换离子为+n 价，则离子交换反应可以用式（1-32）表示，其中R 代表固相上的可交换位点，并设定离子交换平衡常数为Kex，且满足平衡方程式（1-33）。

进一步引入总的离子交换容量Te。假定吸附态的活度为1，则吸附态的活度等于该吸附态所用交换位点数占离子交换容量的当量分数，即式（1-34）；与此同时，各吸附态之间还满足式（1-35）。将固相的交换位点作为组分，利用式（1-35）在质量平衡方程中加入离子交换相，就可以计算固相上的可交换位点被各种离子占据的情况，从而可在任何平衡中加入离子交换平衡。

在实际体系中，往往同时包含多种反应，程序会按照一定的流程设计调用多个反应模型协同解决问题，最终经过多次迭代自检后输出计算结果。

这里就程序中使用的几个模型作简要介绍，以便理解输出结果的含义。表面配位模型将溶质在吸附剂表面上的吸附视为溶液中的离子与吸附剂表面的配位基团间发生的配合反应。吸附模型可以分为经验模型与机理模型。习惯上将传统的吸附模型，如Kd模型、Langmuir模型、Freundlich模型等称为经验模型，而将表面配位模型称作机理模型。经验模型不涉及吸附机理，所得的经验常数只适用于回归这些常数的实验数据所包括的范围，而对于实验条件如p H、Eh、离子强度、被吸附离子及溶液中其他离子浓度的影响难以给出定量的解释，表面配位模型则能从机理上给出较为满意的解释。

吸附剂与溶液接触，吸附剂表面的官能团因解离或吸附溶液中的离子导致其表面带电荷。p H 是影响表面电荷的主要因素，其他阴离子、阳离子的吸附也影响表面的荷电状态。在除零电荷点（point zero charge，PZC）外的其他p H 条件下，表面总是带电荷的。为了保持电中性，在吸附剂表面附近分布着与表面所带电荷符号相反的离子，于是在吸附剂表面便形成了双电层，带电荷的吸附剂表面便具有了一定的电位，我们称之为表面电位。表面附近的离子因受这一电位的作用，其活度与本体溶液中该离子的活度不同。

吸附剂的官能团可用≡SOH 或＞SOH 等类似的符号表示，其中符号“≡”或“＞”表示官能团S—OH 以共价方式与基体物质结合。有的吸附剂的官能团表现为一元酸形式，通过式（1-36）质子化，这类吸附剂称为“1-p K”型吸附剂；有的吸附剂则表现为二元酸形式，其反应式为式（1-37）和式（1-38），这类吸附剂称为“2-p K”型吸附剂。从吸附表面官能团分布来看，有的吸附剂表面官能团分布均匀，吸附性质均一，这类吸附剂称为单位点吸附剂；有的吸附剂表面官能团对于给定离子的吸附能力不均一，有些位点吸附能力较弱，有些位点吸附能力较强，这类吸附剂称为双位点吸附剂。

CHEMSPEC采用的表面配位模型包含非静电模型（NEM）、恒电容模型（CCM）、扩散层模型（DLM）、基本斯特恩模型（BSM）和三层模型（TLM）。不同模型之间的主要区别在于所采用的双电层结构之不同，不同表面配位模型的双电层结构示意图如图1-10所示。

图1-10　不同表面配位模型双电层结构示意图

在非静电模型中，吸附剂表面电荷对表面配位反应没有影响，处于吸附剂表面附近的离子活度与该离子在本体溶液中的活度相等。非静电模型的表面配合反应及相应的平衡常数在形式上与溶液中的配位反应完全相同。

在离子强度比较大（＞0.01 mol/L）的溶液中，扩散双电层可以用一个平行板电容器来近似描述。按照恒电容模型，所有的表面配合物位于平行板电容器的内层（图1-10中的0平面），反号离子则位于外层，背景电解质离子不发生吸附。

在扩散层模型中，吸附剂与被吸附离子的配位反应发生在σ0层，形成内层配合物，背景电解质不发生吸附。扩散层模型认为反号离子同时受到平面电场与热运动的双重作用形成由近到远浓度逐渐减小的扩散层。三层模型是在扩散层模型的基础上，增加了σ1平面，被吸附离子与表面功能团既可生成内层配合物（位于σ0层），也可生成外层配合物（位于σ1层），反号离子位于σd层中。基本斯特恩模型认为在固液相界面存在两个平行的吸附平面，σ0平面位于内部，吸附H+与OH-，被吸附的其他离子位于σ0平面之外的σ1平面，两个平面被一个电容为C1的区域分开，反号离子位于σ1平面之外的扩散层中。

在高放废物地质处置领域，研究的吸附剂多为天然矿物，组成比较复杂，最终的吸附可能是多种模型同时作用的结果。对于复杂吸附体系的研究主要包括广义构成法（GC，generalized composite）与组分加和法（CA，component additivity）。

广义构成法将复杂的吸附视作一种结构单一均匀的整体，将各官能团总的吸附等效为某一特定官能团的吸附结果。一般吸附材料表面往往发生羟基化，将不同的羟基化吸附位点简化为同一种羟基官能团，最后给出总的表面吸附配位常数。这种方法得到的参数比较适合组成固定的体系，但由于其大量简化，需要的信息与实验参数较少，应用得十分广泛。

组分加和法是将吸附剂总的吸附结果看作吸附剂中各单一物质吸附结果经加权后的总和。该方法适合组分相对明确的体系，需要的参数较多，在进行大量模拟后可以探究各单一物质对总吸附量的贡献百分数，预测各物质组成在一定范围内变动时总吸附量的变化。

热力学数据库包含了软件计算运行需要的所有信息。不同的软件都有自己特定的数据格式，数据库的品质对计算结果至关重要。数据库要求具有完备性、准确性、内部一致性、可溯源性、误差满足要求、可读性和易于修改和补充等性质。数据库中应包含体系中所有可能生成的种态。平衡常数和活度系数的可信度会直接影响最终的计算结果。

目前CHEMSPEC采用了具有一定权威性的瑞士PSI的数据库，加上离子交换与表面配位的组分与种态，共有125个组分和555个种态，包含了Th、U、Pu、Np、Am、Tc、I、Se等高放废物地质处置安全评估中关键核素的热力学信息。数据库的扩建通过以下方式进行，收集已经得到实验测定或者经过专家审核的数据，进行类比近似和通过经验公式推导。CHEMSPEC使用的表面配位吸附数据库采用的是德国罗森道夫（Rossendorf）放射化学研究所开发的一个表面吸附热力学专家系统（Rossendorf expert system for surface and sorption thermodynamics，RES3T）中的数据，截至2016年4月30日，该数据库包含了135种矿物，1 668个比表面积测量数据，1 551个表面位点数据记录，4 938个表面配合反应和2 705篇参考文献。总的来看，早期有关矿物的研究较多，目前研究较多的是对早期吸附实验的补充，不同黏土矿上的吸附实验也在逐步开展中。我们需要做的工作是在RES3T系统中根据实际研究体系筛选数据库，按照CHEMSPEC 可以识别的格式输入，就能实现扩展CHEMSPEC 数据库，不断扩大CHEMSPEC 的应用范围。

最近，随着数据库的不断完善、理论模型的改进以及程序计算的优化，CHEMSPEC的计算能力与准确性也有了较大的提升。为使更多的国内同行更好地使用该软件进行有关计算，我们在此重点介绍CHEMSPEC 软件的新进展，并用实例说明软件目前可以计算的问题。

（1）增加了表面配合模型

表面配合模型将吸附剂对溶质的吸附视为吸附剂表面基团与溶液中的溶质离子发生了表面配合反应。表面配合反应与溶液中的配合反应最大的区别是前者涉及固液两相反应，而后者的反应物与生成物均处于液相中。

在表面配合模型中，吸附剂表面的官能团被看作一类组分，称为表面组分，凡组成中包含表面组分的物种称为表面物种。吸附离子与吸附剂吸附位点（表面官能团）相互作用生成的配合物分两种类型：若吸附离子与表面位点以形成化学键的方式结合，形成的表面配合物我们称为内层配合物；若吸附离子与吸附剂吸附位点以静电相互作用形式形成离子对（二者之间常有水分子），这样形成的配合物我们称为外层配合物。为与内层配合物区别，外层配合物通常用连字符“-”或“…”表示，以区分内层与外层，例如

表面配合模型可用多种方法处理静电项，相应地便有多种表面配合模型。CHEMSPEC有5种表面配合模型，分别为：非静电模型、恒电容模型、扩散层模型、三层模型和基本斯特恩模型。具体的处理方法和应用范围可参考CHEMSPEC 用户手册，这里不再赘述。

（2）数据库完善

毋庸置疑，软件计算结果的可信度高度依赖于数据库。一个好的热力学数据库所包含的数据必须具备完备性，准确性和内部一致性。在所给条件下以显著量存在的种态若未包括在相应的数据库中，计算结果将会与实际情况偏离。自开始编写该软件至今，我们对CHEMSPEC的数据库不断进行完善、处理和更新工作。

CHEMSPEC目前收录了PSI/Nagra 和LLNL 两个数据库，用户可以选用。CHEMSPEC还将某些离子交换和表面配合反应的物种包含在内，使得PSI/Nagra数据库包含125个组分和555个物种，LLNL数据库包含265个组分和3 017个物种。这些数据已经转换成CHEMSPEC 的数据格式，即组分文件components.dtb，物种文件species.dtb和化学计量矩阵文件stoichiometry.dtb。由于LLNL的数据库过于庞大，原始数据库也存在诸多问题，因此计算时常常不能收敛。最近，我们花费了非常大的精力重新审查了原始数据库，修正了其中发现的错误，补充了一组关键数据（用于O2分压控制的redox反应计算），并加入了几十组表面配合反应数据。

迄今公开发表的表面配合数据非常多，涉及的吸附剂包括晶态和无定形难溶无机物，天然和合成的高分子材料，天然矿物、土壤及生物物质，研究过的吸附物质包括无机阴、阳离子和有机物。但遗憾的是，对于同一体系，不同作者给出的表面配合反应平衡常数很难进行比较，以至于迄今没有被普遍认同的表面配合反应热力学数据库。造成这种现象的原因是多方面的。①天然矿物的实际组成与结晶状态因产地不同存在差异。②表面配合反应假定溶液/固相界面和固相表面溶液/本体溶液均达到热力学平衡，这点是不能完全保证的，尤其是固液相界面达到热力学平衡可能很慢。在不能达到热力学平衡时，不同作者报告的数据可能来源于不同的时间点。③每种表面配合模型都有一个或多个可调参数。例如不同平面间的电容值C0、C1或C2。不同作者所取的值不同，也会影响表面配合反应的平衡常数。正因如此，德国罗森道夫（Rossendorf）放射化学研究所开发了一个表面和吸附热力学专家系统——RES3T。该专家系统目前包括139个矿物、1 725个比表面积测量数据、1 591个表面位点数据、5 141个表面配合反应和2 842篇参考文献，可以供研究者免费查阅（http：//www.hzdr.de/db/RES3T.login）。

最近一些文献报道了锕系元素和碳酸根离子与碱土金属形成了非常稳定的三元配合物。尤其是钙与碳酸铀酰形成的三元配合物在中性和碱性条件下是铀的优势种态。我们从中选择了一些可信的物种，在数据库中添加了Ba UO2和Ca2UO2（CO3）3等11个三元配合物的有关数据。

向数据文件中手工添加新的数据不但非常麻烦，而且极易出错。为此，我们编写了相应的预处理程序，用户只需将收集到的数据按照规定的格式填写到预处理程序的输入文件中，运行后，新数据就可自动加入相应的数据文件中。

离子和中性分子在溶液中的活度系数计算的准确度对于种态分析计算结果的可靠性影响很大。溶液的离子强度较低时，各种模型差别不大。离子强度＞0.1后，不同活度系数模型给出的计算结果差别明显。对高离子强度溶液的活度计算，Pitzer公式可给出较满意的结果，但需要的参数太多又不易得到。鉴于此，经济合作与发展组织（OECD）推荐使用专属离子相互作用模型（SIT）。我们从文献中收集了尽可能多的SIT参数到CHEMSPEC的数据文件中。

（3）计算优化

使用最新的数据库进行计算时，由于数据库更加完善，搜索到体系可能生成的物种也更加复杂多样，导致软件在计算时十分烦冗，耗时长且不易收敛。为此我们增加了两个物种搜索模式ipick=4和ipick=5。其中，ipick=4是指对搜索到的可能形成的物种，只保留所有的溶液物种，难溶物种全部忽略，但在结果中列出了难溶物种的饱和指数SI。用户根据难溶化合物的SI可以判断在给定的条件下哪些难溶化合物一定不会生成，哪些难溶化合物生成倾向最大，从而为正式计算时对难溶化合物的取舍提供参考。ipick=5是指程序默认接受搜索到的全部溶液物种，将所有的固相物种逐一显示，要求用户对之做出取舍。在进行计算时，可以先令ipick=4，不让固相物种生成，计算速度和收敛性较好，用户从results.dat文件中查到各个固相物种的SI，做选取后再令ipick=5，对固相物种进行取舍。如此，便可在软件计算结果的收敛性和准确性之间做到有效平衡。