Pu在膨润土层中的反应性迁移模拟研究(2)——Pu的种态分布及反应性迁移分析

刘东旭,黄流兴,赵振华,胡立堂,司高华,叶远虑

(1.西北核技术研究所,陕西 西安 710024;2.北京师范大学 水科学研究院,地下水污染控制与修复教育部工程研究中心,北京 100875;3.生态环境部 核与辐射安全中心,北京 100082)

含Pu高放废物的多重屏障处置研究中,需要评估Pu在膨润土等工程屏障材料中的长期迁移影响。基于地下水入渗景象分析,在膨润土回填材料逐渐饱和的过程中,水-岩体系特性可能发生改变和演化,进而影响Pu在膨润土层中的长期迁移。然而,如何采用数值模拟方法获取膨润土演化过程中Pu的迁移规律,仍是国内外地质处置和核素迁移研究的挑战性问题[1]。核素迁移模拟方法研究中的关键是对真实物理-化学过程的理解和模型化,需要在模拟计算中充分考虑水动力迁移与地球化学反应机制,但传统研究方法中,基于对流-弥散方程(CDE)的模型在理论和应用方面均存在局限性,如模型过于简化,仅考虑了水动力弥散作用和基于分配系数Kd的线性平衡吸附机制,而未充分考虑溶解-沉淀、离子交换、表面配位、氧化还原等地球化学反应过程,不能客观揭示反应性元素在地下水系统中的迁移行为[2]。关于Pu等锕系元素的迁移机制,由于氧化态、溶解度效应、胶体机制、矿物表面的吸附和解吸等多种因素影响,目前尚未完全了解特定水-土体系中Pu的种态分布、表面配位模型和迁移转化过程[3]。

反应性迁移模拟(RTM)耦合了地球化学热力学与动力学框架内的流动与迁移过程,它与地球化学模拟不同,后者主要计算静态系统的地球化学平衡、物种形成和热力学状态,并未考虑动力学迁移过程[4-5]。20世纪90年代以来,随着动力学反应理论、实验数据和计算机模拟的发展,TOUGHREACT、OpenGeoSys、CrunchFlow、HPx、CORE2D等RTM软件程序持续发展,其优势是克服了CDE模型过于依赖Kd的不足,真正考虑了对流-弥散-反应过程,这也有助于定量分析工程屏障材料的长期演化、微生物和胶体对核素迁移的影响等关键问题[4-6]。对这些问题的理解和认识将有助于构建更加合理、可靠的性能评价模型,但是目前研究程度还比较有限。RTM研究的挑战性在于:所需数据较复杂且依赖于具体场地和特定材料,研究者往往会面临数据资料不足的情境;反应平衡常数和动力学速率不易获取;核废物处置的环境效应不能通过实验验证,只能尝试模型预测[1,7]。

国际上RTM的应用十分广泛,包括地球关键带低温近地表环境中多相流体流动、溶质运移、地球化学反应等多过程模拟[8],以及CO2地质封存、放射性废物处置等能源与环境问题研究[1,4]。国内在RTM研究方面相对滞后,早期主要是理论方法介绍和技术跟踪,包括水-岩相互作用、耦合水化学模型进行溶质迁移数值模拟等[9-11]。后来,国内学者自主研发的CHEMSPEC[12]以及国际上流行的PHREEQC等程序,用于铀尾矿地下水中U的吸附及反应运移模拟、地下水中砷的水动力-水化学耦合模拟等[13-14]。但在放射性废物处置性能评价研究中,国内更多是基于地球化学模拟来分析不同元素(U、Np、Pu等)在水溶液中的种态分布、吸附和材料演变[15-18],在应用RTM预测分析膨润土工程屏障和围岩材料中Pu的迁移方面研究较少。

由于Pu自身化学行为的复杂性和环境介质的多变性,目前模拟预测Pu的环境迁移行为仍较困难[3]。尽管国内外诸多学者研究了Pu的溶解、配位以及不同矿物(如蒙脱石、黏土、针铁矿等)对Pu的吸附,探讨了pH、Pu氧化态、Fe/Mn矿物等因素对矿物表面上Pu吸附/解吸的综合影响[3,19-21],但极少聚焦于膨润土工程屏障系统中Pu的反应性迁移建模分析。国外仅见少量基于热-水-力-化学(THMC)多场耦合模型分析假想处置情景下Pu在Kunigel-V1膨润土和黏土围岩中迁移的研究[21],该研究采用TOUGHREACT构建膨润土、黏土材料的化学演化模型以及Pu的反应性迁移模型,同时采用表面配位模型来模拟分析Pu的吸附行为,最终得到了处置系统中屏障材料的化学演化规律以及10 000年内Pu的浓度分布。国内司高华[22]利用PHREEQC计算了膨润土工程屏障水环境中Pu的化学形态和溶解度,并在传统CDE模拟框架内基于TOUGH2软件开展了Pu迁移的数值模拟,但并未考虑膨润土化学演化和表面配位吸附机制。

本文结合高放废物处置工程屏障材料性能评估研究,在地下水-膨润土长期演化分析基础上,首先通过地球化学模拟分析处置情景地下水中Pu的种态分布特征;然后,结合文献资料获取Pu在膨润土上吸附的配位反应数据,构建表面配位模型来表征柯尔碱膨润土对Pu的吸附作用;最后,在TOUGHREACT中耦合渗流扩散、溶解-沉淀、表面配位等关键过程,预测分析Pu在柯尔碱膨润土层中的长期迁移行为。通过综合考虑对流-弥散和地球化学反应机制,基于反应性迁移模拟方法进行Pu迁移预测,克服传统CDE模型的局限性,为高放废物处置工程屏障性能评价提供技术依据。

1 材料与方法

1.1 模拟程序及模型设计

综合运用TOUGHREACT、PHREEQC/PhreePlot和Visual MINTEQ,结合动力学模拟和地球化学模拟,建立耦合渗流扩散、溶解-沉淀反应、表面配位吸附、放射性衰变等过程的反应性迁移模型,模拟分析地下水-膨润土地球化学演化过程中Pu的种态分布特征、迁移转化规律。

TOUGHREACT是将地球化学反应引入到TOUGH2而开发的,可用于模拟孔隙-裂隙介质中多相流运动、热运移、溶质运移和水-岩-气化学反应的耦合过程[23]。PHREEQC是美国地质调查局(USGS)开发,用于模拟低温水文地球化学反应过程的计算机程序,具有离子交换、蒸发、流体混合、表面配位吸附、动力学模型、一维迁移和反演模型等模拟功能[24]。PhreePlot是嵌入了PHREEQC程序的化学绘图软件,基于PHREEQC计算结果来生成图像,操作流程相对简单。Visual MINTEQ是一个基于MINTEQA计算内核的化学平衡模型软件,可用于计算自然界水体中的元素形态、氧化-还原、溶解-沉淀、吸附反应等,包含了DLM(扩散层模型)、CCM(恒电容模型)、TLM(三层模型)、CD-MUSIC(电荷分布-多位点表面配位)等表面配位模型[25]。

工程屏障-围岩系统概念模型、TOUGHREACT计算模型、环境数据、边界条件等建模信息参见文献[26]。Pu迁移模拟中,假设在渗流作用下,废物体中核素缓慢释放并扩散迁移,在膨润土层底部设置观测点D,以分析1 m厚膨润土层对Pu的阻滞效果;模型的上、下边界为Dirichlet边界,左、右边界为零流量边界;模型上边界初始压力设定为1.5×105Pa,其余部分初始压力为1.0×105Pa;模型中核素初始浓度为0;239Pu的放射性衰变作用通过其衰变常数来表征。

1.2 表面配位

表面配位模型的一个基本假设是吸附自由能由化学作用和静电作用组成。金属在矿物氧化物表面的吸附主要取决于溶液中带电离子与固-液界面处结合位点之间的静电相互作用,这些与pH相关的结合位点可能处于质子化、中性和去质子化状态。描述矿物表面质子化和去质子化的反应为[27]:

(1)

(2)

上述质子化和去质子化反应分别有一个平衡常数K+和K-,因此又称为2-pK模型。DLM模型固-液界面如图1所示。由图1可见,DLM模型可描述矿物表面羟基与溶液中正负离子发生配位反应以及表面电极附近的净电荷从界面区域到溶液本体的分布。根据Gouy-Chapman理论,DLM模型通过引入考虑静电作用的玻尔兹曼因子来描述平衡常数质量作用方程,如K+可表示为:

(3)

图1 DLM固-液界面示意图

其中:{ }表示物种的活度;F0为法拉第常数,F0=96 485.34 C/mol;ψ0为表面电势,mV;Z为离子的电价;R为气体常数,R=8.31 J/(mol·K);T为绝对温度,K。

DLM考虑了电解质离子在扩散层中的解离,电解质对质子迁移的影响可认为是扩散层电荷与电解质离子作用的结果。常用的表面配位模型都能较好拟合核素的吸附数据,但是DLM扩散层模型所需参数最少,且相关吸附数据较为丰富,而TLM、CD-MUSIC等复杂模型往往存在所需参数较多、适应性较差、不易建模等问题[27]。因此,本文采用2-pK DLM来表征膨润土对核素的吸附作用。

1.3 Pu反应性迁移数据

Pu的化学性质非常复杂,可能以多个氧化态存在。本文通过对吸附过程的综合描述,并结合不同情景设定,间接体现Pu在水溶液中的物种、吸附/解吸、溶解、胶体、金属腐蚀产物等因素对核素迁移的总体影响。

1) Pu在水溶液中的物种分布

通常所考虑的Pu氧化态、水溶液物种及其平衡常数如表1[21,28]所列。这些常数主要取自EQ3/6生成的GEMBOCHS DATACOM.V8.R6数据库,该数据库已用于研究Pu(Ⅳ)和Pu(Ⅴ)在典型膨润土矿物上的表面配位反应[21]。

表1 Pu的水溶液物种及其平衡常数

2) Pu的吸附/解吸

描述吸附/解吸反应的两类常用方法是线性等温吸附模型和表面配位模型。前者采用分配系数Kd来表征吸附作用,进而可简单计算出CDE模型所需的阻滞系数Rd,便于实现建模计算,但也存在过于依赖实验条件、无法合理表征地球化学反应过程等不足;后者是更为严谨的方法,它采用一系列表面反应和表面物种描述固-液相相互作用,能够反映pH、配位反应、不同吸附位点等系统特征,但构建吸附模型所需的表面配位反应及数据相对难获取[1-2,5,7]。文献[21,28]研究了Pu在铝硅酸盐矿物上的吸附及表面配位模型,考察了Pu(Ⅳ)和Pu(Ⅴ)与二氧化硅、三水铝石和典型膨润土/黏土矿物的相互作用,这为详细研究Pu与膨润土矿物相互作用并建立机理模型提供了宝贵的数据资料。参考此研究成果,本文采用表面配位模拟方法来表征Pu的吸附/解吸,利用铝醇和硅醇两类位点(>AlOH和>SiOH)描述Pu在膨润土上的表面配位,这在本质上等同于一般表面配位模型中的弱位点1和2,相关的表面配位反应及平衡常数如表2[21]所列,结合柯尔碱膨润土比表面积实测值与文献[28]数据,得到表面配位模拟中矿物特性参数,如表3所列。

表2 Pu在膨润土表面吸附的表面配位反应及其平衡常数

表3 表面配位模拟中矿物特性参数

3) 固相Pu溶解

通常,固相Pu指的是Pu氧化物/氢氧化物,主要包括Pu(OH)3(s)、Pu(OH)4(am)、PuO2(OH)。科学界虽然对固相Pu的相图了解较少,但对Pu(OH)4(am)或PuO2(hyd,aged)的研究较多,它们表示相同的Pu(Ⅳ)氧化物/氢氧化物[29-30]。研究表明,Pu氢氧化物或胶体会控制Pu的溶解度;除Pu(Ⅳ)的水相沉淀物外,Pu(Ⅳ)水解还会形成胶体,并与Pu(OH)4(am)共同控制Pu溶解[31-32]。然而,在实验中很难区分这些反应及其控制作用,本文暂未考虑Pu胶体的迁移,设定Pu的溶解受Pu(OH)4(am)控制。

2 结果与讨论

2.1 Pu种态分布特征

地下水入渗经过花岗岩、膨润土屏障层后,将与含Pu废物体接触,可能载带Pu并向外迁移。Pu释放迁移模拟时,首先需要考虑Pu在入渗地下水中的氧化态和存在形式。Pu对氧化还原条件敏感,在水环境中存在多种氧化态,同时,Pu具有丰富的配位化学,导致其环境行为较为复杂。Eh-pH图(又称相图或优势区域图)是从酸碱度和电极电势两个方面反映特定元素在环境中的优势形态变化,从该图可获取反应中组分的生成条件以及其稳定存在的范围,广泛应用于地下水中化学反应及平衡问题研究[33]。本文通过地球化学模拟得到Pu在地下水中的Eh-pH图和种态分布,如图2～4所示。

图2 Pu在地下水中的Eh-pH图

图3 不同Eh条件下Pu氧化态的分布(pH=9.5)

图4 不同pH条件下Pu的种态分布

根据文献[26]模拟结果,由于膨润土的pH缓冲性能维持了偏碱性环境,入渗地下水的pH演化范围为8.1～10.3(与膨润土相互作用后pH约为9.5),同时,天然地下水通常呈还原性。由图2～4可见,在上述条件下Pu主要以Pu(Ⅳ)氧化态和Pu(OH)4(aq)形式存在。这与天然水中Pu的种态分布规律相一致[29,34],也反映了研究区特定的地下水、膨润土特征,该结果对于本文工程屏障材料中Pu的迁移评估具有实际意义。

图2～4的模拟结果与文献[35]的研究结果相一致,即pH和Eh对Pu的氧化态、种态形式影响较大,天然环境条件下,Pu几乎不会以离子态存在,且以Pu(Ⅳ)和Pu(Ⅴ)为主,而Pu(Ⅴ)又会被还原为Pu(Ⅳ),使得Pu表观溶解度进一步降低、吸附性增强而迁移性降低。

需要指出,虽然水溶液中Pu(Ⅳ)通常具有难溶、难迁移的特性,然而,Pu(Ⅳ)也易于水解形成胶体或纳米颗粒;在弱碱性条件下,胶体态或聚合态Pu(Ⅳ)将很大程度上影响Pu的氧化还原反应、溶解和迁移性能[30],因此,反应性迁移模拟中将设定不同情景,来模拟分析Pu的长期迁移规律。

2.2 Pu在膨润土层中的吸附阻滞特征

作为工程屏障材料,膨润土对核素起到重要的吸附阻滞作用,因此,如何描述膨润土中Pu的吸附行为至关重要。传统的分配系数Kd方法具有简单、易行、应用广泛的特点,因此,科学界对不同类型膨润土中Pu的Kd进行了广泛研究,结果表明,各种膨润土对Pu普遍具有很强的吸附性,Pu的Kd约为5 800～20 000 mL/g[35-36]。然而,Kd方法存在明显缺陷,它采用恒定的Kd常数来描述复杂的吸附过程,不能反映系统的pH变化、液相配合物的存在形式、吸附相表面不同位点的配位性质等客观因素;相反,表面配位模型考虑了这些影响因素,并可给出更接近实际的等效Kd值[37-38]。

本文基于DLM模型来分析表面配位吸附过程,并得到该模型导出的Kd值,以便于理论模型的推广应用。在考虑Pu(Ⅳ)/Pu(Ⅴ)氧化还原对的情况下,推导得到Pu的Kd与pH关系曲线,如图5所示。由图5可见,pH显著影响Pu的Kd,pH=8时,Kd最大(约16 500 mL/g)。如文献[26],研究区地下水的pH为8.1,该地下水与膨润土发生相互作用后,pH约为9.5,在这种偏碱性的环境条件下,Pu的Kd仍高达约14 500 mL/g,这表明膨润土对Pu具有强烈的吸附作用。通常,实际处置环境中地下水循环缓慢、呈偏碱性和还原性条件,膨润土演化也趋向于维持或增强这种偏碱性的还原环境,因此,可以预计基于表面配位模型的模拟结果将类同于Kd接近最大值、膨润土对Pu吸附阻滞强度很大的情景。Kd接近于0时对应于强酸或强碱性条件,这在天然环境中是不现实的,但本文从研究角度设定一种保守情景,主要考虑Pu在渗流作用下的扩散迁移,假定胶体/纳米颗粒等机制对Pu的迁移具有促进作用,而将表面配位吸附强度减弱1 000倍,通过类比推测可知,此时的等效Kd约为15 mL/g。

图5 Pu的Kd与pH关系曲线

2.3 Pu的反应性迁移规律

废物体中239Pu属于长寿命核素,其半衰期为24 100年,安全评价中通常考虑至少10个半衰期的时间跨度。在如此长的时间尺度下,难以准确预知未来环境条件的变化,核素迁移模拟结果不可避免存在较大不确定性,而多情景预测在一定程度上有助于增强模拟预测的合理性。在上述研究基础上,设定两类情景进行模拟分析:1) 基本情景,即核素迁移过程考虑渗流扩散和表面配位吸附机制;2) 保守情景,即核素迁移过程以渗流扩散为主导而忽略较强的吸附阻滞机制。本文模拟区域中部3 m×3 m子区域为废物体(1 m×1 m)与1 m厚膨润土屏障层,在膨润土层底部D点(x=2.5,z=-4)处观测核素的穿透迁移行为。

基本情景下的模拟结果如图6所示。由图6可见,Pu释放后1年时近似为初始状态,此时Pu浓度晕集中于废物体区域,浓度处于10-9mol/L量级;1万年时,该浓度晕扩散范围略有增大,其中心区域浓度降低约1个量级;24万年时,Pu浓度晕范围又略有增大,但局限于废物体附近的膨润土区域;50万年时,Pu浓度晕衰减至约10-14mol/L量级水平,但仍然没有迁移出1 m厚的膨润土层。由此可见,由于膨润土的低渗透性和强吸附性,Pu的扩散迁移范围很小而将长期滞留于膨润土中。

a——1年;b——1万年;c——24万年;d——50万年

保守情景下的模拟结果如图7所示。由图7可见,Pu浓度晕的变化较明显,源项释放后100年时,该浓度晕已呈现出由中心区域向外径向衰减的特征,Pu在膨润土中的扩散迁移距离约为0.5 m;1万年时,该浓度晕呈持续扩大趋势,大部分区域Pu浓度降低了约1个量级;24万年时,在垂向渗流作用下,该浓度晕在径向扩散衰减的同时略呈下垂状,且扩散迁移范围进一步扩大出膨润土层;50万年时,由于没有考虑吸附机制,在渗流作用下,Pu浓度晕呈垂向运移趋势,其中心区域的峰值浓度已扩散至膨润土层底部。由此可见,假如没有膨润土的强吸附滞留作用,即使在低渗透条件下,Pu的扩散迁移范围仍然较大,1 m厚的膨润土层已不足以实现含Pu废物的长期有效隔离;这也说明,如果存在这一保守情景(如在胶体、微生物等作用下,Pu与膨润土之间的表面吸附作用失效,Pu的迁移受控于渗流和扩散作用),那么,Pu的长期迁移存在风险性。正因如此,当前科学界已广泛重视可能导致渗流扩散主导的迁移机制,如处置情景下各种膨润土、黏土材料的侵蚀以及所产生胶体的迁移性、胶体促进核素迁移的理论模拟等研究[36,39-42]。

a——100年;b——1万年;c——24万年;d——50万年

基本情景和保守情景下,膨润土层底部D点处Pu的液相浓度变化如图8所示。由图8可见,两种情景下,膨润土层底部D点处Pu的穿透浓度差别很大。50万年尺度下,基本情景中D点Pu浓度很小,可认为Pu不会穿透膨润土层。保守情景中Pu浓度短期内持续增大至10-9mol/L量级水平,直至数万年后,该浓度才呈衰减下降趋势,这进一步表明,由于膨润土的低渗透性和强吸附性,Pu的扩散迁移范围非常有限,1 m厚的膨润土层能满足安全性要求,而一旦失去该吸附滞留作用,膨润土层将失去屏障效能,因此,实际中仍需重视渗流扩散主导的反应性迁移。

图8 基本情景和保守情景下Pu在膨润土层底部D点处的液相浓度变化

3 结论

1) 结合柯尔碱膨润土矿物组分、地下水化学特征等数据资料,通过地球化学模拟得到了Pu的Eh-pH图、氧化态及种态分布特征,分析认为,研究条件下Pu主要以Pu(Ⅳ)氧化态和Pu(OH)4(aq)形式存在。

2) 采用DLM扩散层模型和文献数据建立了Pu在柯尔碱膨润土中吸附的表面配位模型,导出了Pu的Kd-pH关系曲线;在pH=8或9.5的基本情景下,等效Kd高达16 500 mL/g或14 500 mL/g,表明膨润土对Pu的吸附阻滞作用很强。

3) Pu长期迁移模拟预测表明:基本情景中,由于柯尔碱膨润土的低渗透性和强吸附性,即使经历24万年乃至50万年,Pu浓度晕仅局限于废物体附近的回填屏障层区域,1 m厚的膨润土层可以实现阻滞Pu迁移的安全性目标;保守情景下,由于忽略了强吸附作用,Pu在垂向渗流作用下呈向下扩散迁移趋势,最终在相对短的时期内可穿透该膨润土层,使得膨润土失去屏障作用。鉴于此,建议重视膨润土胶体行为、Pu反应性迁移机制等研究,并考虑Pu一旦穿越出膨润土屏障层的后果分析,以增强性能评价的科学性。

感谢北京师范大学王金生教授、吉林大学许天福教授和田海龙教授在数值模拟方面给予的支持。