放射性核素在地下水中迁移问题一直以来都是国内外研究的热点问题。铀尾矿库作为整个铀矿冶设施主要放射性污染源项，其包容了大量天然长寿命放射性核素和非放射性有害物。在雨水淋滤作用下，这些物质会逐渐向地下水环境中发生迁移，通过地下水运移、弥散进入到人类生存环境，最终通过饮水途径对人类形成放射性内照射。因此，了解和掌握在役和退役铀尾矿库内放射性核素在地下水中迁移现状和迁移规律，对铀矿山环境影响评估和矿区环境污染的整治具有重要意义。

1 核素迁移模拟研究现状

目前，地下水中核素迁移研究主要集中于室内实验室研究、地下实验室研究、天然类比研究、计算机模拟研究等方面^[1]。由于铀尾矿库设计安全年限一般较长(IAEA为1000年)，长时间范围内难以通过常规实验手段进行核素时空变化规律分析，因此，尾矿库内放射性核素在地下水中迁移研究多通过计算机数值模拟实现，数值模拟方法正逐渐成为定量评价地下水污染的有效工具。目前地下水中核素迁移模拟研究中，理论较为成熟、应用广泛的污染物输运模型主要有水动力学模型和多组分水化学平衡模型。其中，水动力学模型是将溶质在地下水中运移看成对流、弥散、衰变、生物降解等因素综合作用的结果，侧重研究水动力场下核素时空分布规律。多组分水化学平衡模型，又称水文地球化学模型，侧重研究溶质在地下水中水化学作用下组分、形态分布、迁移规律。对于铀尾矿库内放射性核素在地下水中运移数值模拟，目前多采用美国地质局开发MODFLOW软件，基于流体水动力学模型、有限差分数值解法，用于孔隙介质、裂隙介质中的地下水流动和核素迁移问题^[2]。

2 数学计算模型以及求解方法 2.1 地下水水动力学模型

地下水系统中任何污染物质的运移、预测，都是以地下水动力学为基础，了解放射性核素与土壤介质的吸附-解吸作用、放射性核素与地下水之间相互作用因素，分析核素在地下水系统中流体动力弥散机理，是进行放射性核素迁移数值模拟基础。在地下水中溶质输运过程中，影响放射性核素迁移主要有4个因素。

2.1.1 对流输运

对流输运(又称平流输运)，是指地下水在含水地层介质孔隙中流动时，携带着水中污染物以相同运移速度在孔隙中迁移。

2.1.2 弥散作用

地下水中污染物弥散由机械弥散和分子扩散两部分组成。其中，机械弥散主要是由于流体通过多孔介质时，速度不均匀而造成物质输运；分子扩散则是由于溶质浓度分布不均匀，溶质在浓度梯度作用下由高浓度向低浓度处运动，使其浓度分布均匀，是分子布朗运动的一种表现^[3]。

2.1.3 吸附作用

污染物在介质颗粒表面上会发生吸附从而阻滞其迁移，使其在多孔介质中的输运速度小于地下水的运移速度，使得污染迁移范围减小。吸附滞留包括吸附解吸、离子交换、可逆性沉淀、过滤等廓清作用。

2.1.4 衰变降解作用

放射性核素在地下水中输运时，受核素原子核内部运动规律影响，核素自发发生放射性衰变，其衰变服从指数分布规律。

考虑稳定均匀流场条件下、对流、弥散、吸附、衰变降解作用基理，依据地下水水动力学模型，可建立放射性核素三维迁移的基本方程：

$ \begin{align} &\frac{\partial C}{\partial t}={{K}_{x}}\frac{{{\partial }^{2}}C}{\partial {{x}^{2}}}+{{K}_{y}}\frac{{{\partial }^{2}}C}{\partial {{y}^{2}}}+{{K}_{z}}\frac{{{\partial }^{2}}C}{\partial {{z}^{2}}}-{{\mu }_{x}}\frac{\partial C}{\partial x}-{{\mu }_{y}}\frac{\partial C}{\partial y} \\ &\ \ \ \ \ \ \ \ \ -{{\mu }_{z}}\frac{\partial C}{\partial z}-\frac{{{\rho }_{b}}}{{{\eta }_{e}}}\frac{\partial \left( {{K}_{d}}C \right)}{\partial t}-\lambda C~ \\ \end{align} $

1)

其中，C-地下水中放射性核素污染物浓度；K_x、K_y、K_z-分别为地下水流动方向、横截方向、铅垂方向上的弥散系数；μ-地下水水流速度；ρ_b-多孔介质密度；η_e-有效孔隙度；K_d-分配系数。放射性核素三维迁移方程描述了地下水动力场中的任一点任一时刻的水动力参数、介质参数与浓度之间的关系，方程中各项表征物理含义分别为：等号左边项表示放射性核素浓度随时间的变化，等号右边前3项表征水动力弥散作用造成的核素迁移；等式右边4~6项为对流项，等式右边最后2项目分别表征吸附作用项、衰变降解项。

2.2 有限差分网格数值解法

由于放射性核素三维迁移方程为抛物线型二阶线性偏微分方程，只有在理想条件下，才能求得解析解；一般情况下，放射性核素三维迁移方程很难求得解析解，只能通过数值方法求得近似解。本次模拟采用的MODFLOW软件，采用有限差分网格法求解放射性核素三维迁移方程，其基本原理为：将时间和空间区域划分成若干网格，用末知函数值在网格结点上的值所构成的差分近似代替所用偏微分方程中出现的各阶偏导，从而把表示变量连续变化关系的偏微分方程离散为有限个代数方程，然后解此线性代数方程组，以求出放射性核素在各网格结点上不同时刻的浓度^[4]。

3 实例研究 3.1 研究区概况

研究区铀尾矿库位于我国南方丘陵地区，是一个有着30多年运行历史的铀尾矿砂处置场所，库内现贮存尾砂1.9×10³万t^[5]。铀尾矿库近似呈南宽北窄的梯形，尾矿坝体由九个拦水坝段和丘陵山头组成。其中，尾矿库南北方向长约1.4 km，东西方向宽为0.9~1.5 km。尾矿库北面、南面、西面被某河流环绕，河流距离库区最近距离约5 km，常年平均河水位为48 m。库区海拔高度为62~104 m，地形高差30~40 m，地貌为河流Ⅲ级阶地，由上自下主要发育为亚粘土层、卵石层、粉砂岩层及砂岩层，库区地下水潜水含水层受大气降水影响，直接进行补给，向河流进行排泄。尾矿库所在区域出露地层主要为第三系古新统东塘组、始新统霞流市组茶山坳段、第四系更新统和全线统，主要岩性为红层碎屑岩类及第四系松散沉积物。区内属亚热带潮湿气候，冬寒期短，炎热期长，雨量充沛。常年平均气温为17.5℃~17.9℃，年均降水量1233~1363 mm，年均蒸发量1394~1533 mm，相对湿度77%~79%。

3.2 模型的概化与建立

根据铀尾矿库周边区域水文地质条件，本次取铀尾矿库区所在区域分水岭到河流区段作为研究区域。由于铀尾矿库区地表高程、地下水水位标高均高于周边区域上游地区、下游地区，尾矿库3面被河流所环绕，模型概化时选取河流作为河流边界。铀尾矿库东侧2~3 km为地下水分水岭，将其设置隔水边界。区内地下水总体流动方向为自西向东方向，因此可将研究区的南北两侧设置隔水边界。受大气降水因素影响，将整个研究区设置降水入渗边界。在尾矿库北侧边缘有一水文监测井。为了便于建立数学模型的建立与求解，本次铀矿库区地下水中流动可假定研究区地下水流动为二维稳定流、孔隙介质为均匀多孔介质、放射性核素在多孔介质中的吸附-解析过程为可逆过程，基于放射性核素多孔介质水动力学溶质运移模型原理，考虑放射性核素在地下水中对流、弥散、岩土介质对核素吸附、核素衰变等作用的综合影响，建立铀尾矿库地下水中放射性核素方程。

放射性核素迁移方程：

$ \begin{align} &{{R}_{d}}\frac{\partial C}{\partial t}={{K}_{x}}\frac{{{\partial }^{2}}C}{\partial {{x}^{2}}}+{{K}_{y}}\frac{{{\partial }^{2}}C}{\partial {{y}^{2}}}-{{\mu }_{x}}\frac{\partial C}{\partial x}-{{\mu }_{y}}\frac{\partial C}{\partial y} \\ &\ \ \ \ \ \ -\lambda C{{R}_{d}}+{{C}_{0}}Q{{\eta }_{e}}~ \\ \end{align} $

2)

初始条件：

$ 铀矿库区:C\left( x, y, 0 \right){{|}_{\left( x, y \right)\in {{\mathit{\Omega }}_{0}}}}={{C}_{0}}, 0\le t\le {{t}_{1}}~ $

3)

$ 铀尾矿库区以外其他区域:C\left( x, y, 0 \right){{|}_{\left( x, y \right)\notin {{\mathit{\Omega} }_{0}}}}=0~ $

4)

边界条件：

$ 河流边界: H\left( x, y, t \right){{|}_{{{\mathit{\Gamma }}_{0}}}}={{H}_{0}}~ $

5)

$ 地下水分水岭隔水边界 H\left( x, y, t \right){{|}_{{{\mathit{\Gamma }}_{1}}}}={{H}_{1}} $

6)

其中，C-地下水中放射性核素污染物浓度，Bq/L；K_x、K_y-分别为地下水纵向、横截方向弥散系数，m²/s；μ_x、μ_y-地下水纵向、横截方向水流速度，m/s；λ-核素衰变常数，s^-1；R_d -放射性核素阻滞系数，无量纲；R_d -放射性核素初始浓度，Bq/L；Q-单位时间向单位体积含水层抽注水量，s^-1；H₀、H₁-地下水位高程，m。

本次实例研究选取铀尾矿库在内的11 000×8 800 m区域为研究范围，模拟的放射性核素²³⁸U(Ⅵ)、²²⁶Ra在地下水中二维迁移。利用MODFLOW软件建立模型时，依据研究区地质以及水文地质资料，将铀尾矿库研究区域地层从上到下依次分为亚粘土层(含粘土层)、卵石层、粉砂岩层，各层水动力学参数见表 1^[6]。每层采集204个高程点，然后进行克力格插值，生成各层层面。为了便于求解放射性核素迁移方程，将地层每层划分为350×300网格，每个网格代表距离为31.4 m×29.3 m，再利用有限差分公式进行求解，将上述方程变成线性方程组，通过强隐式方法迭代进行求得近似解。根据铀尾矿库提供的监测资料，铀尾矿库运营期间废水中²³⁸U浓度、²²⁶Ra分别为8.75、1.74 Bq/L，其放射性衰变常数分别为4.92×10^-18 s^-1、1.37×10^-11 s^-1。研究区降雨入渗系数取0.058，年降雨入渗量取75 mm。河流、地下水分水岭隔水边界地下水高程分别取48、100 m。

3.3 模拟结果与分析

本次选取²³⁸U、²²⁶Ra核素作为模拟对象，相关参数见表 2。依据地下水中放射性核素迁移方程，在MODFLOW软件中建立起相关数学模型，对模型进行概化，输入相关计算参数，采用差分网格计算方法，研究铀尾矿库放射性核素²³⁸U、²²⁶Ra活度浓度在时间、空间尺度上的变化。数值模拟时，依据铀尾矿库场地及周边环境地下水流场，分别模拟时间为30、50、100、200、500、1000 a条件下铀尾矿库区²³⁸U、²²⁶Ra活度浓度。监测井中²³⁸U、²²⁶Ra活度浓度与时间变化关系见图 1、图 2。

由图 1、图 2模拟结果可知，该铀尾矿库退役后0~100 a期间，尾矿库中²³⁸U、²²⁶Ra活度浓度呈快速减少趋势。铀尾矿库退役100 a后，铀尾矿库北侧边缘监测井中²³⁸U、²²⁶Ra活度浓度分别为7.55、0.95 Bq/L，其活度浓度分别降低了14%、45%。这主要是由于早期铀尾矿库内含有大量中低放射性的尾液，其含水量较高，尾液中放射性核素经过包气带土层阻滞、吸附等作用后，随尾液或降水进入到地下水潜水含水层，在地下水动力场作用下，随地下水流动向周边环境迁移，其迁移速度较快，导致监测井中地下水中²³⁸U、²²⁶Ra活度浓度快速减少；100~500 a期间，铀尾矿库中尾矿砂中含水率逐渐降低，导致尾矿库中放射性核素迁移速度放缓，监测井中²³⁸U、²²⁶Ra活度浓度逐渐减少。铀尾矿库退役500 a后监测井中²³⁸U、²²⁶Ra活度浓度分别为6.18、0.76 Bq/L，与初始浓度相比，²³⁸U活度浓度降低了29%，²²⁶Ra活度浓度降低56%；500~1000 a期间，铀尾矿库内²³⁸U、²²⁶Ra迁移速度极小，监测井中²³⁸U、²²⁶Ra活度浓度基本上无变化。模拟结果表明，随着时间的向后推移，铀尾矿库中²³⁸U、²²⁶Ra核素会逐渐向周边环境迁移，其活度浓度呈现先快速减少、后趋于稳定的变化趋势。

IAEA规定铀尾矿库退役后应至少满足1000年的安全要求，本次对研究区铀尾矿库模拟了1000年后的地下水中²³⁸U、²²⁶Ra活度浓度分布，模拟结果分别见图 3、图 4。

由图 3、图 4可知，研究区尾矿库退役1000 a后，铀尾矿库地下水中放射性核素²³⁸U向西迁移约0.5 km，向北迁移约1.7 km，向南迁移1 km，放射性核素²³⁸U迁移面积约达到6 km²；地下水中放射性核素²²⁶Ra向西迁移约0.5 km，向北迁移约0.7 km，向南迁移0.7 km，放射性核素²²⁶Ra迁移范围为3.8 km²。模拟结果表明，铀矿库退役1000 a后，尾矿库中放射性核素²³⁸U、²²⁶Ra在地下水流场作用、核素的弥散作用、沿途多孔介质对核素²³⁸U的吸附阻滞作用以及核素本身的衰变作用下，在靠近铀尾矿库侧核素²³⁸U、²²⁶Ra活度浓度相对较高，随着相对距离的增加，核素²³⁸U、²²⁶Ra活度浓度逐渐减少。

4 结论与建议

依据研究区铀尾矿库的水文、地质、气象等条件，基于地下水动力学模型，将放射性核素在地下水中的迁移概化二维条件下的稳定流动，利用MODFLOW软件，借助有限差分网格计算方法，模拟得到放射性核素随时间、空间的变化特征。模拟结果表明：研究区尾矿库内放射性核素²³⁸U、²²⁶Ra活度浓度随时间逐渐减少，1000 a后其浓度达到一稳定值，放射性核素²³⁸U、²²⁶Ra在地下水中逐渐发生迁移，导致库区1~2 km距离内放射性活度浓度升高，可能对库区周边地下水环境产生影响。核素在地下水介质中迁移受多种因素影响，本次数值模拟仅简单考虑地下水动力场作用、放射性核素本身衰变因素，未考虑在地下水中运移过程中地下水化学组分、PH值、化学形态、腐殖层、胶体等因素影响，因此建议下一步加强地下水地球化学方面的研究，更好地应用于核素迁移模拟。