氡示踪河底地下水-河水相互交换数值模拟研究

＊乐庆国

（香港城市大学 香港 999077）

1.概念模型

(1)河水-地下水交换模型

天然河流里，具有周期的三角沙丘河床地形相当普遍。以河流流向为X轴，垂直河床面为Y轴，建立二维河床剖面模型。模型中，拟定河床沉积物中存在连续的三角形沙丘微地形，缓坡方向与河流流向相同，陡坡方向与河流流向相反。边界条件中，以模型左边界为速度入口，模型右边界为压力出口，河床、河面为无滑移边界，河床由5个连续的三角形沙丘构成，每个沙丘波长1m，位于0.75m处有一尖峰，高为0.115m。采用ANSYS Fluent模拟河流紊流流动，采用的剖分方法为三角剖分，节点数在90000个左右，尽可能避免由于网格剖分数量给模拟结果带来影响，以此模型模拟计算河水给其下伏的潜流带的压力及其分布。

(2)河流紊流模型

ANSYS Fluent模拟河流紊流原理为有限体积法求解不可压缩流体纳维-斯托克斯方程（Navier-Stocks Equation）。所用到的方程为RANS方程和k-ω湍流模型，RANS方程如下：

式中，ρ、μ为流体的密度和动力黏度；Ui（i=1,2）、ui（i=1,2）分别为流体的平均速度和瞬时速度；g为重力加速度；P为流体平均压强；Sij为应变速率张量，其式为：

RANS方程与k湍流动能和湍流耗散率ω有关，其式为：

其中，τij为雷诺切应力；vt为涡流运动黏度；δij为克罗内克符号。

湍流动能κ方程为：

湍流耗散率ω方程为：

其中，μt=ρvt，α=5/9，β=0.075，β*=0.09，σκ=σω=0.5。

(3)水流和溶质运移模型

河床形态引起的潜流交换主要包括泵吸交换（Pumping Exchange）和冲淤交换（Turnover Exchange）。假设三角沙丘河床不发生移动，潜流交换主要由泵吸作用引起。孔隙介质饱和渗流可用方程描述：

假设河床沉积物为均质各向同性，h为水头；K为渗透系数；n为孔隙度。

222Rn在孔隙介质中运移涉及对流和弥散，且由于其具有放射性，它的产生和衰变方程为：

2.河流流速对氡示踪潜流交换的影响

(1)模拟方案

水动力条件作为潜流带内物质和能量交换的重要驱动力，不仅控制着地表水和地下水的水量交换，影响两者的混合比例，同时控制着随水流进入潜流带内各物质的停留时间、混合运移和时空分布模式。在本研究中，河床面初始压力与河流流速直接相关联，不同的河流流速会产生不同的河床面初始压力，初始压力在COMSOL Multiphysics中以压力边界的形式存在，压力边界是边界条件的一种，是确定一个水文地质模型能否求解的一个关键参数。本研究选定河流流速为变量，定量研究河流流速对初始压力的影响，并探究后者对于222Rn示踪潜流交换的影响。模拟方案的取值见表1，其余参数恒定不变，见表2。

表1 各个模型河流流速取值

表2 其余参数取值

(2)结果与分析

图1 不同河流流速下河床的压力分布

①河流流速变化下河床压力分布

从图2来看，四个模型都展现出了河水在一定流速驱动下，由于沙丘形态，流速横截面缩小，冲击着迎水面的沙丘，在大约沙丘迎水面的中部压力最大，河水由此进入河床沉积物，在沙丘的左端和右端压力最小，河水在此返回河流。沙丘迎水面和背水面的压力差随着流速的增加而呈平方增长每一个沙丘的压力线型十分相似，但是在总体趋势上呈现出压力随着河水流动距离增加而减少，河水每流经一个沙丘，沙丘压力峰值下降大约20%，最小值则下降大约50%。值得注意的是，在第三个沙丘的右侧0.75m处低凹区，四个模型都出现了负压。此外，模型1，2在沙丘顶0.75m处都出现了显著低于左右侧的压力值，模型3，4则不明显或基本可以忽略。河床面的压力值与河流流速的平方呈正相关关系，模型2和模型4的河水流速只相差一个量级，但河床压力数值却相差两个量级，这意味着河流流速增加导致的河床面压力增加的边际效应十分显著。

图2 不同河流流速下潜流带的水头分布图

由于不同河流流速下每个沙丘的压力线型相对一致，四个模型的水头分布图较为相似，具体表现为：每个沙丘中部0.5m处水头较大，作为系统的源；沙丘左右侧0～0.1m和0.9～1m水头较小，作为系统的汇；较浅的区域水头差较大，随着沉积物深度增加，水力梯度逐渐减小，潜流交换强度随着沉积物深度增加而降低。由于模型3，4的第三个沙丘的低压力区相较于前两个沙丘压力差过于巨大，低压力区已蔓延到模型底部，意味着来自前两个沙丘的河水更有可能会深入潜流带底部与深部地下水进行交换，再由模型右侧的低压力区离开。水头数值与初始压力呈正相关关系。

②河流流速的变化对氡示踪潜流交换的影响

图4至图6可以看出，其相同点是，在模拟的早期，222Rn示踪效果较好，往往在模拟时间t=1d时，河水就已扩散到潜流带中相当大的一部分，这部分区域可以用低浓度222Rn示踪出来，即该区域的222Rn浓度远低于地下水中的222Rn浓度，约1至2个量级；越往深处，由于含水层固体颗粒矿物（226Ra，222Rn的母体）衰变产生222Rn进入潜流带，潜流带的222Rn浓度逐渐升高，逐渐接近地下水氡浓度，氡示踪潜流的交换的效果逐渐减弱。

模拟期第10天已基本展现出稳定的态势，相较于第20天模拟期结束的最后结果相差不大，甚至完全一致；由上到下，从低浓度区到高浓度区，每一个模型都出现了一层渐变过渡带，近似呈“3V”型，高浓度区的氡经由这一渐变带扩散至低浓度区，并在此过程中浓度逐步降低。原因为222Rn经历了数个半衰期，从模拟一开始时进入模型中的222Rn在抵达浓度渐变带区域基本衰变殆尽，对于浓度渐变带及其之下的深度，222Rn对于示踪潜流交换的效果不佳。

针对不同区域，从图3到图6，可以清楚的发现初始压力和氡浓度分布之间的关系：初始压力越大，低浓度的222Rn区域就越大。模型1对应极低流速下的河流，例如，河流下游接近入海口的河段。由于水动力条件较为薄弱，在渗透系数不太大的情况下，河水向潜流带扩散有限，仅能影响河床以下约0.1m区域，浓度渐变带约0.3m，对于222Rn来说，示踪再深处的区域则力不能及。模型2展现的氡浓度分布图相较模型1，低222Rn浓度区域明显大得多，大约0.2m宽，渐变带宽度约0.4m。模型3，4结果则不相同，在模拟期结束的第20天，低222Rn浓度区域都扩散到了近模型底部，只是模型4的扩散速度更快，它们的低浓度区域和浓度渐变带也较模型1，2更厚，约0.5m。

图3 模型1河流流速为0.06m/s 222Rn浓度分布图

图4 模型2河流流速为0.2m/s 222Rn浓度分布图

图5 模型3河流流速为0.5m/s 222Rn浓度分布图

图6 模型4河流流速为1m/s 222Rn浓度分布图

造成模拟结果如此分异的原因为河床面上的压力差不同，模型1，2的沙丘迎水面—背水面压力差较小，难以驱使河水深入潜流带进行交换；而模型3，4的压力差高达数百帕，这有助于提高潜流交换强度，以对流为主的溶质迁移速率也会加强，在浓度图上的低浓度区域也相对较大。而在这种情况下，222Rn就能够很好地示踪模型中的潜流交换。

3.结论

基于COMSOL Multiphysics建立的河水与地下水水流-溶质运移耦合模型模拟结果，表明在模拟的早期，潜流带立即形成一片由河水构成的低浓度区，随着河水进一步向深处扩散，氡的运移减慢，低浓度区扩张逐渐放缓稳定，并从上至下形成一条浓度由低变高的渐变带。ANSYS Fluent模拟河流紊流的结果表明，地表水流条件和微地形形态控制着初始压力，其中初始压力对于前者尤其敏感，与前者的平方呈正相关；后者则会改变初始压力的分布。二者皆对氡示踪潜流交换具有显著影响，河水水流增大，河床地形坡度增大，氡示踪效果相对更加显著。

水文地质参数对于氡示踪潜流交换的效果具有重要影响。渗透系数K≤0.1，介质渗透性较差，模型中低浓度区域仅局限于水沙界面，222Rn的示踪效果不佳；如果K≥1，低浓度区域显著大于前者，且渗透系数越大，低浓度区域越大，222Rn的示踪效果越为明显。孔隙度影响着多孔介质的放射性元素衰变产生的氡进入地下水的速率，孔隙度越大则氡的产生速率越大，越大的孔隙度会使氡的浓度平衡深度变浅。弥散系数控制着在进入含水层的河水在对流路径上的弥散，弥散系数越大，河水中的氡能迁移的更远、更深，低浓度区域也越大，222Rn的示踪效果也会越好。地下水氡浓度值对于氡的示踪影响甚微，其值越小，越难以与地表水区别开，222Rn的示踪效果也越差。通过对以上地表水条件和水文地质参数变化影响的研究，刻画了每个参数对于氡示踪地表水—地下水的交换的影响，可为地下水与地表水的交叉研究和实践提供理论基础。