放射性核素迁移的包气带-饱和带耦合模拟研究

张 雪，胡继华，施晓文，王小元

(中国核电工程有限公司，北京 100840)

为保证核设施建设、运行及退役等阶段土壤及地下水环境的安全，根据国家核安全法规、标准和导则的要求，必须通过实验和数值模拟等手段，充分掌握核设施厂址地区放射性核素在土壤及水环境中的迁移扩散规律，提供工程设计所需的参数资料。

以往对假想事故条件下放射性液体的地下水环境影响评价，假设放射性液体瞬时进入饱和地下水环境，不考虑包气带对核素迁移的迟滞作用。对于放射性核素来说，包气带对放射性核素的阻滞作用非常明显，李书绅等[1]的研究表明，在自然条件下，可吸附核素的迁移距离非常短，因此研究放射性核素在包气带中迁移是必要的。放射性核素在包气带中的迁移主要以垂直向为主，且与地下水环评的评价范围相比，水平向的运动几乎可以忽略。

将包气带水分运动处理为一维垂向运动，饱和带地下水为三维运动，将两种不同维数运动的问题进行耦合而得到拟三维数值模型。包气带一维垂向模拟选用Hydrus-1D软件，饱和带三维模拟选用GMS软件。将潜水面作为包气带模型下边界，同时作为饱和带模型上边界，通过包气带模拟获得包气带下边界地下水流量及放射性核素浓度随时间变化值，并将其作为饱和带的源汇项(上边界条件)输入到饱和带模型，实现包气带-饱和带的松散耦合。

1 国内外研究进展

大多数饱和-非饱和水流耦合模型都基于MODFLOW建立，基于一维非饱和水流运动和三维地下水运动的含水层数值模拟系统如图1所示。

图1 基于MODFLOW建立的拟三维饱和-非饱和水流系统

Havard等人[2]建立了LINKFLOW模型，该模型将模拟非饱和区垂向流动的Richards方程的有限差分解与模拟饱和区水流的三维有限差分模型MODFLOW进行耦合。Facchi等人(2004)[3]将土壤水运动模型SVAT与MODFLOW进行了耦合。Niswonger和Prudic (2004)[4]、Niswonger等人(2006)[5]在MODFLOW-2005 的基础上发展了非饱和水流运动计算程序包UZF1。Twarakavi等人(2008)[6]采用类似的方法建立了一个更为完善的饱和-非饱和水流运动模型，采用HYDRUS-1D程序包模拟非饱和带水流运动，并与MODFLOW-2000进行耦合。Yakirevich等人(1998)[7]采用有限差分法建立了拟三维饱和-非饱和水流运动和溶质运移耦合模型。Kuznetsov等人(2012)[8]进一步发展了Yakirevich的耦合模型，采用完全三维水流运动方程描述地下水与毛细管区水流运动,并与非饱和一维水流运动方程进行耦合。

林琳等(2005)[9]根据多孔介质中区域饱和-非饱和水分运动的特征，建立基于迦辽金有限单元法的拟三维模型。朱炎(2013)[10]将非饱和带水流运动简化为一维垂向运动，保持饱和带三维运动的特点，将两者进行耦合。查元源(2014)[11]将Richards方程模型及水均衡模型与三维地下水模型耦合，建立了区域饱和-非饱和三维水分运动耦合模型。于洋等(2017)[12]通过GMS和Hydrus的结合，实现了包气带-饱和带水盐运移的耦合模拟。张雅楠(2019)[13]明晰了HYDRUS-MODFLOW模型的耦合原理及运行机制。

2 研究区概况

研究区属大陆性干旱气候，降水少，蒸发强烈。地势北高南低，北部为基岩山地，海拔1 250.0～1 440.0 m，南部为碎石戈壁，海拔一般为1 229.0～1 275.0 m。多年平均降雨量62.9 mm，多年平均蒸发量为2 338.9 mm。地下水可分为松散岩类孔隙水和基岩裂隙水两类。松散岩类孔隙水区域渗透系数为2.78×10-4～7.20×10-3cm/s；基岩裂隙水区域渗透系数为3.97×10-6～1.53×10-4cm/s。地下水水位埋深在2.80 m～20.06 m之间，水位标高在1 205.25 m～1 227.8 m之间。

3 包气带模拟研究

3.1 数值模型的建立

包气带一维垂向模拟选用Hydrus-1D软件，建立模型时做如下假设：包气带中水流和核素为垂向一维流，核素的释放源在地表；包气带介质对核素的吸附是可逆的，为一级线性反应；只存在一级化学反应，且不考虑核素的衰变产物；忽略温度和气体的影响。

试验场地包气带厚度约3 m，假设罐体中的放射性液体发生泄漏，并由地表垂直入渗补给包气带，泄漏总量为8.50 m3，核素种类包括C-14、Se-79、Pd-107、Cs-135和Sn-126五种，活度浓度分别为9.40×104Bq/L、8.80×103Bq/L、1.50×104Bq/L、5.60×104Bq/L和6.64×104Bq/L。

3.2 初始条件

(1)假设初始条件下水流处于稳定状态，水流模型初始值由压力水头给出：h(z,0)=h0(z)，-H≤z≤0。式中：h为土壤水压力冰头，H为包气带厚度，z为垂直坐标，设潜水面处压力水头为0，向上递减，地表处压力水头为-3 m。

(2)核素迁移模型初始时刻选在放射性液体刚从地表泄漏，包气带中放射性液体初始浓度设为0。

3.3 边界条件

(1)水流模型上边界为变水头/流量边界:

(1)

泄漏会产生地表积水，t1为泄漏持续时间，H(t)为地表积水厚度。即地表有放射性液体泄漏时(0≤t≤t1)，上边界为随时间变化的水头边界；泄漏停止后(t>t1)上边界为0流量边界。k(h)为非饱和水力传导度。

下边界选择潜水面，设为定水头边界：

h(z,t)=0，z=-H，t≥0

(2)

(2)核素迁移模型上边界为定浓度通量边界：

(3)

式中，θ为体积含水率，D为水动力弥散系数，v为水流垂向流速，c为核素在包气带中的浓度，qs(t)为地表水分通量，C为泄漏液体中放射性核素浓度。

下边界为零浓度梯度边界。

3.4 模型主要参数确定

土-水特征曲线参数θs、θr、n、α和Ks通过现场土-水特征曲线试验和现场包气带渗水试验获取初值，并利用包气带水分迁移试验，通过模型识别验证获取模型参数取值；干密度ρd取测试值；弥散度aL选择多孔示踪试验值。参数取值列于表1。

表1 参数取值表

表2 分配系数和衰变常数取值

3.5 模拟结果及分析

根据包气带-饱和带耦合模拟思路，包气带的模拟结果将作为饱和带的源汇项(上边界条件)输入到饱和带模型中，模拟得到潜水面处水流速度及C-14、Se-79、Pd-107、Cs-135和Sn-126五种核素浓度随时间变化如图2～图7所示。

图2 潜水面处水流速度随时间变化

图3 潜水面处C-14活度浓度随时间变化

图4 潜水面处Se-79活度浓度随时间变化

图5 潜水面处Pd-107活度浓度随时间变化

图6 潜水面处Cs-135活度浓度随时间变化

图7 潜水面处Sn-126活度浓度随时间变化

由图2～7模拟数据结果可知：

(1)包气带对水流运动(核素迁移)存在延迟作用：放射性液体并非瞬时进入饱和带，包气带具有一定的渗透能力，水流只能以一定的速率渗入饱和地下水环境(本次模拟完成水流入渗大约需要30 d)，核素主要在对流作用的驱动下渗入饱和地下水环境，也即包气带对核素迁移存在延迟作用。

(2)包气带对核素迁移存在滞留作用：由于包气带介质对核素的吸附作用，使迁移至潜水面处核素的活度浓度缓慢增大而并非瞬时增大。对于吸附性较强的核素(Se-79、Pd-107、Cs-135和Sn-126)，其迁移至潜水面处的核素活度浓度远小于泄漏的放射性液体核素活度浓度，Se-79、Pd-107、Cs-135和Sn-126在潜水面处的活度浓度峰值分别为3.67 Bq/L、705 Bq/L、3.67 Bq/L和0.045 Bq/L，对应的年待积有效剂量分别为8.19×10-4mSv、0.018 9 mSv、8.19×10-4mSv和5.83×10-5mSv，已小于WHO规定的0.1 mSv。而对于吸附性较弱的C-14，迁移至潜水面处的核素活度浓度与泄漏液核素活度浓度差别不大，C-14在潜水面处的峰值活度浓度为9.37×104Bq/L，为泄漏液活度浓度的99.68%，对应的年待积有效剂量为15.50 mSv，与泄漏的放射性液体(15.55 mSv)相比只降低了0.05 mSv。

即包气带对核素的迁移存在明显的迟滞作用，根据包气带模拟结果，接下来只对潜水面处核素活度浓度对应年待积有效剂量大于评价准则(0.1 mSv)的C-14开展耦合模拟研究。

4 饱和带地下水水流数值模型

4.1 模型范围、边界条件与水力特征

模拟范围与边界条件充分考虑评价区水文地质条件确定。

模拟区内地下水运动符合达西定律。地下水系统的垂向运动主要是大气降水在地表入渗后向下补给含水层，孔隙水含水层垂向补给基岩裂隙水，模拟区内没有明显的隔水层，根据地下水位长期监测结果，该地区地下水水位随季节变化不明显。因此，模拟区地下水系统可概化为非均质三维稳态地下水流系统。

4.2 源汇项处理及水文地质参数设置

模拟区地下水主要接受大气降水和田间灌溉水补给，主要排泄项为蒸发和开采。

(1)降雨入渗：厂址区年均降雨量为62.9 mm，降水入渗补给系数根据地表岩性参考水文地质手册推荐的经验值。

(2)蒸发：由地下水位观测资料，水位埋深大于极限蒸发深度4 m，故不考虑饱和地下水的蒸发。

(3)灌溉和开采：模拟区内180眼农业灌溉井的地下水年开采量约4 360万m3，计算得单井开采量约664 m3/d，单井灌溉面积约199 800 m2，灌溉入渗系数参考水文地质手册的推荐值，取0.2，灌溉补给速率为(664 m3/d×0.2)÷199 800 m2=6.64×10-4m/d。即开采量取664(m3/d)/眼，灌溉补给速率取6.64×10-4m/d。

(4)渗透系数：渗透系数根据水文地质调查及水文地质试验成果，结合含水层岩性，同时参考经验值给出，并在模型识别验证中进行微调。

4.3 模型识别和验证

建立地下水流模型后，通过拟合-校正进行模型识别，本次模型识别与校正标准包括以下三个方面。

4.3.1识别后的模型概化和参数应符合实际的水文地质条件

校核出的含水层渗透系数和降雨入渗补给量分区如图8和图9所示。综合考虑以往研究成果和现场试验数据及经验值，各水文地质参数基本符合实际特征。

图8 渗透系数取值分区图

图9 降雨入渗补给量分区图

4.3.2稳定流状态下观测井地下水水位与计算水位基本一致

本次模拟拟合9眼长期观测孔，观测地下水水位取一个水文年水位平均值，经模型识别后的观测点水位计算值和实测值拟合结果见图10。图中显示实测水位与计算水位基本一致。

图10 观测孔水位拟合图

4.3.3模拟地下水流场与调查得到的模拟区地下水流场基本一致

将模拟的地下水流场和实测的地下水等值线进行对比，结果示于图11。图11中黑色线为根据实测水位数据绘制的地下水等值线，彩色线表示模拟水位分布。从对比图可知，地下水模拟流场与实测流场形态基本一致，地下水流向相同。

图11 地下水等水位线对比图

5 包气带-饱和带耦合模拟

5.1 核素迁移概念模型

核素迁移概念模型可概化为放射性液体通过包气带垂直进入饱和地下水环境后，在对流、弥散、吸附和衰变作用下，随地下水向环境敏感点迁移，核素迁移概念模型如图12所示。

图12 核素迁移概念模型示意图

5.2 包气带一维垂向迁移模型与饱和带三维模型的耦合

将包气带模拟获得的潜水面处水流速度及C-14浓度随时间变化值作为上边界条件输入到GMS饱和带模型。

5.3 包气带-饱和带耦合模拟结果及分析

C-14模拟结果示于图13。由图13可知C-14在对流作用下不断向下游迁移，在弥散作用影响下，其污染范围逐渐增大，在地下水稀释及衰变、吸附等作用下，C-14活度浓度逐渐降低。事故发生后第1、5和10年，C-14放射性液体前缘迁移距离分别为1.5 m、4 m和4.2 m。

活度浓度单位为Bq/L，图中网格非模型网格，间距1 m。

根据模拟数据结果，C-14在事故工况下发生泄漏后经过包气带渗入饱和地下水环境，在地下水稀释及吸附、衰变、对流、弥散等综合作用下，其活度浓度大大降低，地下水中C-14的峰值活度浓度为2.15 Bq/L，对应的年食入待积有效剂量为9.1×10-4mSv，远小于剂量限值0.1 mSv。事故工况后，放射性液体C-14活度浓度时空变化规律列于表3。

表3 事故工况后放射性液体C-14活度浓度时空变化规律

6 不考虑包气带的模拟研究

为了进一步说明包气带对核素迁移的迟滞作用，开展不考虑包气带的模拟研究(以下简称饱水带模拟研究)，假设放射性液体瞬时进入饱和地下水环境，并将模拟结果与包气带-饱和带耦合模拟结果进行对比分析。

6.1 模拟条件

饱水带模型与耦合模型的区别在于放射性液体进入饱和地下水中的方式：耦合模型中放射性液体经过包气带缓慢渗入饱和地下水，饱水带模型中放射性液体瞬时进入饱和地下水，亦即核素迁移模型的初始和源汇项(上边界条件)设置。

6.1.1核素迁移模型初始条件

耦合模型中C-14初始活度浓度为0 Bq/L。而饱水带模型中C-14初始活度浓度为9.4×104Bq/L，体积为8.50 m3。

6.1.2核素迁移模型源汇项

耦合模型中C-14补给活度浓度随时间变化，详见表3。而饱水带模型将C-14活度浓度作为初始条件，不再考虑其补给源项。

6.2 模拟结果及分析

模拟结果示于图14。由图14和模拟数据结果可知：包气带对放射性核素的迁移存在明显的迟滞作用，在包气带的影响下，地下水中C-14的活度浓度大大降低，放射性液体的迁移距离也明显变小。

活度浓度单位为Bq/L，图中网格非模型网格，间距10 m。

(1)在考虑包气带时，地下水中C-14峰值活度浓度为2.15 Bq/L，对应的公众食入年待积有效剂量为9.1×10-4mSv，小于评价准则0.1 mSv；在不考虑包气带时，地下水中C-14峰值活度浓度为9.4×104Bq/L，对应的年食入待积有效剂量为39.8 mSv，远大于0.1 mSv。

(2)在考虑包气带时，事故发生后第1、5和10年，C-14放射性液体迁移距离分别约1.5 m、4 m和4.2 m；在不考虑包气带时，事故发生后第1、5和10年，C-14放射性液体迁移距离分别约3 m、12 m和17 m。详见表4。

表4 包气带对C-14活度浓度和迁移距离的影响

7 结论

为了能够更加准确地开展地下水环境影响评价工作，本文采用拟三维数值方法开展包气带-饱和带迁移耦合模拟模型研究工作：包气带水分运动处理为一维垂向运动，饱和带地下水为三维运动。为了进一步说明包气带对核素迁移的迟滞作用，开展了不考虑包气带的模拟研究，并将模拟结果与包气带-饱和带耦合模拟结果进行对比分析。

主要结论如下：

(1)包气带模拟研究表明，由于包气带对核素迁移的迟滞作用(延迟和滞留)，对于吸附性较强的核素(Se-79、Pd-107、Cs-135和Sn-126)，其迁移至潜水面处的核素活度浓度远小于泄漏的放射性液体核素活度浓度，对应的公众食入年待积有效剂量已小于剂量限值0.1 mSv，只有吸附性较弱的C-14在潜水面处的核素活度浓度对应年食入待积有效剂量大于0.1 mSv，需对其开展耦合模拟研究工作。

(2)耦合模拟研究表明，C-14在对流作用下不断向下游迁移，在弥散作用影响下，其污染范围逐渐增大，在地下水稀释及衰变、吸附等作用下，C-14活度浓度逐渐降低。地下水中C-14的峰值活度浓度为2.15 Bq/L，对应的年食入待积有效剂量为9.1×10-4mSv，远小于剂量限值0.1 mSv。事故发生后第1、2和10年，C-14放射性液体前缘迁移距离分别为1.5 m、4 m和4.2 m。

(3)饱水带模拟结果表明，包气带对放射性核素的迁移存在明显的迟滞作用，在包气带的影响下，地下水中C-14的活度浓度大大降低，放射性液体的迁移距离也明显变小：在考虑包气带时，地下水中C-14峰值活度浓度对应的公众食入年待积有效剂量为9.1×10-4mSv，小于评价准则0.1 mSv，而不考虑包气带时，地下水中C-14峰值活度浓度对应的年食入待积有效剂量为39.8 mSv，远大于0.1 mSv；在考虑包气带时，事故发生后第1、5和10年，C-14放射性液体迁移距离分别约1.5 m、4 m和4.2 m，而不考虑包气带时，事故发生后第1、5和10年，C-14放射性液体迁移距离分别约3 m、12 m和17 m。