在退役的铀尾矿库中，尾矿砂受到大气降雨的淋滤和地表径流的冲洗溶解进入浅层地下水，其中，核素铀主要溶解态U(VI)的存在对人类健康和周围生态环境造成了严重的威胁(刘奇等，2007).含水层中矿物的吸附影响着U(VI)在浅层地下水中的迁移，对于铀尾矿库的退役治理和环境保护具有重要的参考意义.

吸附对U(VI)的阻滞作用一直以来都备受关注.早期吸附作用的研究一般采用线性吸附等温线模型，马腾等(2000)计算出分配系数后应用MT3DMS(Zheng et al., 1999)对铀尾矿地区的U(VI)迁移进行了研究，该方法将复杂的化学系统简化为仅由单组分吸附离子在固、液相中的分配系数(K_d)来控制的经验过程(Krupka et al., 1999).但实际上由于溶质与围岩矿物的动态反应，水化学组分总是随时间和空间不断发生改变，而且各组分之间的反应是相互影响的(易树平等，2011).其它经验吸附等温线模型有Freundlich和Langmuir模型等，尽管它们在实际中得到了较广泛的运用，但都不能够从机制上反映出其它组分的变化对吸附的影响.

将吸附从化学反应机制上模拟是通过表面络合模型(Surface Complexation Model，SCM)来实现的.Dzombak等(1990)关于水合铁氧化物(Hydrous Ferric Oxide)的吸附作用给出了系统的表面络合模型介绍.Kohler等(1996)经过一系列U(VI)吸附柱实验后提出了U(VI)-石英矿物表面络合反应系列模型.在U(VI)-石英系统中Kohler把石英表面概化为3个类型：只有一种吸附位的石英表面；有强、弱两种吸附位的石英表面反应体系；有强、中、弱3种吸附位的石英表面反应体系.第2类和第3类石英表面又分别各自概化出3种不同的表面反应体系，加上第1类石英表面反应共得到7种不同的表面络合反应模型.其中，第4种双吸附位模型与柱实验运移数据的拟合最好，即第2类双吸附位(强、弱)三表面络合反应式模型，并在此模型基础上推测了pH变化对U(VI)迁移的影响.

为模拟含水层U(VI)迁移中化学组分之间的相互作用，刘奇等(2007)应用水化学模拟软件PHREEQC结合表面络合模型在简单一维水流条件下对U(VI)的吸附作用进行了研究，模拟了溶液-介面之间的复杂化学系统，但PHREEQC无法模拟复杂地下水流条件下的反应溶质运移.基于此，本研究中使用的反应运移软件PHT3D处理了复杂水流条件下的吸附反应迁移问题，并对几种不同的吸附模型进行了比较.

饱和含水层中地下水流数学方程(薛禹群等，2010)为：

式中，K_ii为坐标轴x_i方向与各向异性介质主方向平行时的渗透系数值(m · d^-1)；h为水头(m)；q_s为每单位体积中的源/汇流量(d^-1)；S_s为贮水率(m^-1)；t为时间(d).

溶质组分运移数学方程(Prommer et al., 2010)中，对于液相迁移组分有：

对不可流动组分(如矿物相)：

式中，v_i为沿x_i方向的孔隙水实际平均流速(m · d^-1)；D_ii为流速方向与弥散张量主方向平行时的水动力弥散系数值(m² · d^-1)；r_reac，n为由化学反应引起的浓度变化(mol · L^-1)；θ 为孔隙度，无因次；C^s_n为源汇项中第n种组分的总水相浓度(当为汇项时C^s_n=C_n)(mol · L^-1)；C_n为第n种组分的总水相浓度(mol · L^-1).

在化学反应体系中并非所有组分都要进行详细的描述和模拟(Dai et al., 2004)，通常选择N_m种元素或物质作为主要组分Q^m_i(master species)，而次要组分Q^s_j(secondary species)可以表示为主要组分的线性组合(Yeh et al., 1991)：

式中，v_ji为组成第j种次要组分的第i种主要组分的配比系数；N_s为次要组分的数目.

各组分之间的平衡反应如水相络合、沉淀溶解、氧化还原及离子交换吸附等均满足质量作用方程：

式中，a代表组分的活度，等于组分浓度与活度系数γ的乘积，活度系数可由Davies方程或Debye-Hückle方程算得；K_j为生成第j种次要组分的离子缔合反应平衡常数(Parkhurst et al., 1999).不同类型的反应应进一步推导其质量作用方程具体表达式.

这样在化学反应体系中未知变量就只有各主要组分浓度，因此，根据各主要组分浓度的监测值，再加上水的活度、质量守恒与电荷平衡，就可以得出整个化学反应体系的状态.

线性吸附等温线假设固体基质吸附的溶质的量C^*(mol · kg^-1)与溶质的浓度C(mol · L^-1)之比为常数K_d，即分配系数(L · kg^-1).Freundlich吸附等温线是一种更为普遍的吸附等温线，其非线性关系式为：

式中，K、N为常数.Freundlich吸附等温线与线性等温线均假定固体基质的吸附容量无限大，因而吸附浓度可以随溶质浓度无限增大(郑春苗等，2009).

Langmuir吸附等温线建立在固体表面的吸附点有限的基础上，当所有的吸附点都充满后，固体表面就不再从溶液中吸附溶质(Fetter，2011)，Langmuir吸附等温线形式为：

式中，α为与结合能有关的吸附常数(L · mol^-1)；β为固体所能吸附的溶质的最大量(mol · kg^-1).由此可做出C/C^*与C变化曲线拟合得到两个常数.当溶液浓度低时，Langmuir等温线与线性等温线近似；当溶液浓度高时，吸附浓度达到最大β.因此，Langmuir等温线比Freundlich等温线更能真实地反映高浓度时的野外情形(郑春苗等，2009).

图 1(Fig. 1)

图 1 Freundlich等温线(a)与Langmuir等温线(b)示意图 Fig. 1 Freundlich isotherm(a) and Langmuir isotherm(b)

在表面络合模型中，金属阳离子M²⁺的吸附反应(马腾等，2000)可写为：

表面络合模型以双电层理论为基础，将固体表面看作一种聚合酸，其表面羟基可以发生表面络合反应，并用与处理溶液中络合反应类似的质量作用方程来处理表面过程，但在表面络合平衡过程中需要考虑邻近基团电荷的影响(Goldberg et al., 2007).在质量作用方程上比液相络合作用多一个静电项e^－zFΦ/RT，表面络合反应的表观常数(刘奇等，2007)为：

式中，K^int为表面电荷为零时固有的表面常数；e^－zFΦ/RT称为Boltzmann因子，无因次；Φ为吸附表面电位(V)；z为离子电荷(C)；F为法拉第常数；R为理想气体常数，取值为8.314 J · mol^-1 · K^-1；T为绝对温度(K).

相对于经验模型而言，表面络合模型考虑了电荷对吸附离子和吸附表面的影响，能够直观地解释液相组分的变化对吸附反应的影响(Davis et al., 1978).它包含了一系列不同化学特征和可变参数的模型，这些模型之间的主要区别在于固-液界面的结构特征，如吸附离子的位置及用来描述表界面静电的电荷-电势关系(Waite et al., 1994).

PHT3D由Prommer等(2010)开发，它继承了水流和溶质运移软件MT3DMS及地球化学反应软件PHREEQC(Parkhurst et al., 1999)的功能，其程序模块示意图见图 2.PHT3D可处理复杂(二维、三维)水动力系统下的水岩相互作用中复杂的平衡反应，如水相络合、矿物溶解、离子交换和表面络合等，也可模拟不可移动组分的动力学反应，尤其是微生物活动和NAPL相的溶解(Appelo et al., 2010).

图 2(Fig. 2)

图 2 PHT3D模拟示意图 Fig. 2 PHT3D conceptual sketch

PHT3D在对流弥散和化学反应方程组求解中通过算子分离方法分为两个序列进行，将化学反应项从对流、弥散和源/汇项中独立出来，并在每个时刻上不相互迭代，即序列非迭代法(SNIA)(郑春苗等，2009).反应运移方程的时间差分格式可写为：

式中，C_n^k与C_n^k+1为第n种组分在k与k+1时刻的总水相浓度，L(C_n)^k+1/2为k+1/2时刻的对流、弥散微分算子，Δt为时间步长.在求解多组分的连锁反应时分两步：第一步，由对流、弥散模块(MT3DMS)求解仅有物理运移的k+1时刻浓度C_n^transp：

第二步，紧随第一步完成后按照网格单元顺序独立求解化学反应项，最后得到C_n^k+1：

常用的反应运移模拟软件还有低温序列水文地球化学模拟软件PHREEQC、多相组分平衡模拟软件MINTEQ2(Allison et al., 1991)、多相流溶质运移及适用于高温高压条件下的反应运移软件TOUGHREACT(Xu et al., 1999)，以及变饱和水流及水化学模拟软件HYDRUS-HP2(Šim[AKu。D]nek et al., 2012).

研究场地为一废弃的铀尾矿库，其最早由Yeh等(1991)简化作为研究对象，后来Šim[AKu。D]nek等(2012)做了进一步改进和完善，以使其与实际更加接近.图 3为库区附近地形的二维垂直剖面图，虚线方框内为铀尾矿砂存放范围，长约16 m，地面以下深约8 m.铀尾矿库靠近地表，向右地形逐渐变为斜坡，斜坡底部为一条河流，河流附近地下水位常年基本保持在4 m左右.当地大气降水产生的地表径流下渗到尾矿库和周围土壤中，Yeh给出的尾矿库中渗漏速率为0.0139 m · d^-1，地表降雨入渗速率为0.00139 m · d^-1.尾矿库区的地下水位在12~13 m之间，尾矿库区的海拔高于右侧斜坡下游的河流，而且降雨入渗速率大，潜水流向河流一侧，地下水位逐渐降低.

图 3(Fig. 3)

图 3 研究场地剖面图及地下水等水位线图 Fig. 3 Conceptual model profile and groundwater contour

降水入渗使得库区尾矿砂不断得到淋滤，尽管后来矿砂被转移，但其中溶解态的U(VI)渗漏到地下水中成为污染源.当含有U(VI)的地下水到达河流时，也会引起河流的污染.根据前人大量研究，U(VI)在迁移过程中会受到含水层中某些矿物的吸附，如石英矿物(Kohler et al., 1996)、水合铁氧化物(Dzombak et al., 1990)等.矿物吸附使U(VI)的运移受到阻滞，相反解吸时会促进其运移.

研究场地含水层概化为均匀介质(Yeh et al., 1991)，孔隙度为0.3，水平方向渗透系数为0.378 m · d^-1，水平与垂向渗透系数比值为3，纵向及横向弥散度分别为2.5、0.25 m.

根据水文地质概念模型，选择计算区大小为长105 m，高24 m.网格剖分为24单元×6单元，左侧垂直边界及底部水平边界设为隔水边界，下游河流为定水头边界，上部为降水入渗补给边界.模拟时间为800 d，时间步长为0.5 d.

模拟化学反应体系定义了35种水相物质和14种矿物(Šim[AKu。D]nek et al., 2012)，含水层中定义主要组分为Ca²⁺、CO^2-₃、UO²⁺₂、PO^3-₄、SO^2-₄、Fe²⁺和H⁺.含水层中发生的化学反应包括液相络合、沉淀溶解、吸附解吸等，表 1中给出了化学反应体系中一部分U(VI)的液相络合反应式及平衡常数.含水层初始地下水组分和尾矿库内各组分浓度及补给水流各组分浓度如表 2所示.方解石矿物和石膏矿物保持溶解平衡，矿区外方解石初始含量为4.7×10^-4 mol · L^-1，石膏含量为0；矿区内方解石含量为0，石膏初始含量为3.7×10^-3 mol · L^-1.其它12种次生矿物包括菱铁矿、氢氧化钙、铁铀云母及柱铀矿等，它们的矿物相初始含量为0 mol · L^-1，并处于沉淀溶解平衡状态.

表 1(Table 1)

表 1 化学反应体系-液相络合及表面络合反应(Yeh et al., 1991；Kohler et al., 1996) Table 1 Aqueous complexation and surface complexation formulations in the reaction network

表 1 化学反应体系-液相络合及表面络合反应(Yeh et al., 1991；Kohler et al., 1996) Table 1 Aqueous complexation and surface complexation formulations in the reaction network 化学反应体系 平衡常数logK 表面吸附点位浓度/(mol·L-1) U(VI)液相 UO2+2 + H2O = H+ +(UO2)(OH)+ -5.3 络合反应 2 UO2+2 + 2 H2O = 2 H+ +(UO2)2(OH)2+2 -5.68 3 UO2+2 + 4 H2O = 4 H+ +(UO2)3(OH)42+ -11.88 3 UO2+2 + 5 H2O = 5 H+ +(UO2)3(OH)+5 -15.82 4 UO2+2 + 7 H2O = 7 H+ +(UO2)4(OH)+7 -21.9 3 UO2+2 + 7 H2O = 7 H+ +(UO2)3(OH)-7 -28.34 UO2+2 + CO2-3=(UO2)(CO3)(aq) 9.65 UO2+2 + 2 CO2-3 =(UO2)(CO3)2-2 17.08 UO2+2 + 3 CO2-3 =(UO2)(CO3)4-3 21.7 2 UO2+2 + CO2-3 + 3 H2O = 3 H+ +(UO2)2(CO3)(OH)-3 -1.18 UO2+2 + SO2-4 =(UO2)(SO4)(aq) 2.95 UO2+2 + 2 SO2-4=(UO2)(SO4)2-2 4.0 2 H+ + UO2+2 + PO3-4 = H2(UO2)(PO4)+ 23.2 3 H+ + UO2+2 + PO3-4 = H3(UO2)(PO4)2+ 22.9 Ca2+ + 4 H+ + UO2+2 + 2 PO3-4 = CaH4(UO2)(PO4)2+2 45.24 Ca2+ + 5 H++ UO2+2 + 2 PO3-4 = CaH5(UO2)(PO4)3+2 46.0 表面络合反应 S_aOH + UO2+2 + H2O = S_aOUO2OH + 2 H+ -4.9503 5.197×10-4 S_bOH + UO2+2 + H2O = S_bOUO2OH + 2 H+ -3.4703 1.256×10-5 S_bOH + UO2+2 = S_bOUO2+ + H+ 0.7766 1.256×10-5 S_cOH + UO2+2 + H2O = S_cOUO2OH + 2 H+ -1.0603 1.092×10-7

表 2(Table 2)

表 2 模拟区初始及补给组分浓度(Yeh et al., 1991) Table 2 The initial and recharge component concentration

表 2 模拟区初始及补给组分浓度(Yeh et al., 1991) Table 2 The initial and recharge component concentration 组分 (价态) 组分浓度/(mol·L-1) 含水层 尾矿库 补给 Ca(2) 1.0×10-2 1.0×10-2 1.0×10-3 C(4) 1.5×10-3 1.0×10-2 4.57×10-3 U(6) 1.0×10-7 5.0×10-4 1.0×10-8 P(5) 1.0×10-6 1.0×10-6 1.0×10-6 S(6) 1.036×10-2 2.786×10-1 1.0×10-4 Fe(2) 1.0×10-3 2.0×10-1 1.0×10-3 H(1) 1.0×10-7 3.5×10-2 1.0×10-7

研究中表面吸附模型选自于Kohler等(1996)提出的U(VI)-石英矿物表面络合反应模型系列中的第5种C5模型，含有强、中、弱3种吸附位，该模型具有较长的拖尾效应，吸附解吸效果显著.表 1给出了C5模型3种吸附位的表面络合反应式及表观平衡常数.

经验吸附模拟根据表面络合模拟结果得到观测点的C^*-C曲线，拟合相关参数得到前述提到的几种经验吸附模型，再将其应用到反应运移数值模拟中，比较不同吸附模型的模拟结果.

图 4为不考虑吸附(左)U(VI)反应的运移结果与使用表面络合模型模拟吸附(右)U(VI)运移结果.当t=100 d时，不考虑吸附的U(VI)运移5×10^-7 mol · L^-1浓度线和5×10^-6 mol · L^-1浓度线均超过水平方向40 m，而表面络合中都没有超过40 m.当t=300 d时，不考虑吸附的U(VI)运移5×10^-7 mol · L^-1浓度线已到达水平方向83 m，5×10^-6 mol · L^-1浓度线也达到67 m，而表面络合中两条浓度线仅到达60 m左右.当t=600 d时，两者的运移情况大致相似.前述结果说明了表面吸附对U(VI)的运移有明显阻滞作用，增加了放射性核素铀在地下含水层中的停留时间，使得污染治理难度加大，威胁人类健康和周围生态环境.

图 4(Fig. 4)

图 4 不考虑吸附与表面络合吸附的U(VI)运移浓度变化等值线图 Fig. 4 U(VI)transport without adsorption and with SCM considered

图 5对表面络合模拟结果从分配系数K_d角度进行了分析，并结合图 6给出的U(VI)吸附过程中的pH变化，可以看出反映吸附能力的分配系数K_d与pH有着密切关系：U(VI)运移300 d相对于100 d的pH值在水平方向上逐渐降低，同时分配系数K_d也有相同的变化；并且K_d发生较大范围变化的区域基本与pH发生较大变化的区域重合.这与马腾等(2000)的研究相符合，根据马腾等的研究，在碱性条件下，UO²⁺₂与CO^2-₃、HCO^-₃可以形成不易被表面吸附的碳酸盐络合物，吸附变化受液相浓度变化的影响很小，而在中性或偏酸性条件下吸附变化程度才较大.

图 5(Fig. 5)

图 5 表面络合中的Kd变化 Fig. 5 Kd variation in SCM simulation

图 6(Fig. 6)

图 6 表面络合中的pH变化 Fig. 6 pH variation in SCM simulation

由表面络合中K_d变化可知，在强酸条件下吸附程度变化也不大(图 5、图 7)，观测点A在模拟期基本处于强酸条件下，其吸附U(VI)的量很小(图 8).观测点B、C的pH变化较大，其吸附变化程度也较大.这也解释了图 4中t=600 d时吸附与非吸附模拟结果相似的原因，首先U(VI)浓度受水流和弥散的影响变小，其次较低的pH环境下吸附量就更少，此时的吸附作用可忽略不计.由图 5及图 7三个观测点的K_d变化可以明显看出，K_d在空间和时间上是变化的，而这种变化正是由于反应运移过程中不同地球化学条件下的吸附能力不同引起的.

图 7(Fig. 7)

图 7 3个观测点的表面络合Kd变化 Fig. 7 Kd variation on 3 observations in SCM simulation

图 8为应用表面络合模型进行模拟3个观测点的U(VI)吸附浓度C^*与溶液浓度C的变化曲线，可知吸附浓度并非随溶液浓度的增加而增加或保持最大，相反超过某一浓度值后会减小，这是由于高浓度U(VI)所处的强酸性地球化学条件使得固体基质的吸附能力下降.为使用线性吸附模型，对图 8左边有直线上升趋势的观测点进行拟合，得到分配系数K_d为0.1 L · kg^-1.图 9为C/C^*与C变化曲线，使用Langmuir吸附模型拟合，得到最大吸附容量β为1.0×10^-6 mol · kg^-1，吸附常数α选为7.0×10⁴ L · mol^-1.

图 8(Fig. 8)

图 8 表面络合模拟中U(VI)吸附浓度C*与液相浓度C变化 Fig. 8 Adsorbed C* and aqueous C in SCM simulation

图 9(Fig. 9)

图 9 表面络合模拟中C/C*与C变化 Fig. 9 Variation of C/C* and C in SCM simulation

图 10为使用上述线性吸附模型和Langmuir吸附模型的模拟结果.当t=100 d时，应用两种经验模型得到U(VI)浓度大于3.0×10^-4 mol · L^-1的面积均大于表面络合模型，这是由于线性吸附等温线和Langmuir吸附等温线在高浓度时的高吸附能力，并在很短的时间达到平衡，U(VI)吸附量达到最大，当溶液浓度扩散降低时固相吸附的U(VI)会不断释放出来从而在污染源处能够一直保持较高的浓度；而表面络合模型在高浓度U(VI)时的吸附能力并不是最大(图 5、图 8观测点A)，这是由于尾矿区的低pH强酸性条件抑制了表面吸附反应的进行，U(VI)浓度随对流弥散作用降低迅速，与不考虑吸附的情形大致相当.

图 10(Fig. 10)

图 10 应用线性吸附与Langmuir吸附等温线模型模拟U(VI)运移运移浓度变化等值线图 Fig. 10 U(VI)transport with linear isotherm and Langmuir isotherm

当t=300 d时，线性吸附模拟U(VI)的1.0×10^-4 mol · L^-1浓度锋面到达水平方向35 m位置，Langmuir吸附模拟的浓度锋面到达38 m，均滞后于表面络合模拟的45 m.当t=600 d时，应用两种经验吸附模型的结果为U(VI)浓度大于1.0×10^-4 mol · L^-1的面积(红色虚线内)比表面络合模型的大，位置也相对滞后，其中心在40 m处，高度约为10 m.表面络合模型结果U(VI)浓度大于1.0×10^-4 mol · L^-1的区域中心位置在55 m处，高度小于5 m.经验吸附模型前期U(VI)高吸附浓度导致后期解吸浓度也相对较高，位置也相对滞后.

从上述比较结果来看，线性吸附等温线模型与Langmuir吸附等温线模型的缺点就在于其经验性：两者的吸附能力都仅由U(VI)液相浓度决定，没有考虑到体系中其它组分对吸附的影响.而表面络合模拟结果表明这种由溶液中其它组分对吸附的影响是不可忽视的，不考虑这种影响有时甚至会得出与实际相反的结果.表面络合模型从吸附的反应机制出发建立起来的，其吸附能力能够与周围的地球化学条件相符合，并随地球化学条件的变化而变化.

首先对考虑吸附和不考虑吸附的地下水中U(VI)反应运移进行模拟，发现两者结果差异明显，吸附对U(VI)运移的阻滞作用是明显的，忽略吸附不能够反映真实的场地情况.然后比较了使用不同吸附模型的结果，表面络合模型、线性吸附等温线模型及Langmuir吸附等温线模型三者模拟结果均有差别，表面络合模型能够适应不同的水文地球化学环境，而后两者经验模型其吸附能力仅与单相组分U(VI)的浓度有关，没有考虑到其它组分对吸附的影响.

模拟中的不足之处在于水流模型都是独立于运移模型而不受影响，实际上尤其对于反应运移中如沉淀溶解作用会影响到介质的渗透性，使水动力条件发生变化，此时水流模型要重新建立.