含水层非均质性对潜流带氮迁移转化的影响机制

凡倩莹,刘 洋,冯永林,周继强,王全荣*

(1.中国地质大学(武汉)环境学院,湖北 武汉 430078;2.山东省地质科学研究院,山东 济南 250013;3.自然资源部高寒干旱区矿山地质环境修复工程技术创新中心,甘肃 兰州 730000)

潜流带是氮削减和河流水体净化的主要场所,其水文地球化学过程是影响地下水环境演化过程的主控因素。近年来,工农业活动产生大量的氮污染物[1],导致河流氮负荷逐年增加,河流水体富营养化[2-3]。河流下方及附近沉积物中河水与地下水的混合区(潜流带)可以有效地降低或去除水环境中过量的氮[4-6]。自然水文或人为事件如洪水、降雨、潮汐作用、大坝释水等引起的河水位瞬时改变[7],又进一步增强了河水与地下水的混合,加剧了潜流带溶质的生物地球化学反应。

高频率的河水位波动会促进河流与潜流带间水量[8]、溶质的交换[9]和迁移转化[10-14],使得物质浓度和反应速率[15]在短距离和时间尺度内发生变化,进而影响硝酸盐的去除[15-17]。数值模拟作为水文地质学研究领域中常用的研究手段,因具有有效、快捷、灵活、相对廉价性等优点,已被诸多学者用于研究河水位波动对潜流带生物地球化学过程的影响。Gu等[16]首先采用数值模拟手段研究了河水位波动对潜流带氮循环的影响,但研究中忽略了沉积物中原位碳释放过程产生的溶解有机碳,也没有考虑硝化作用;Shuai等[17]在Gu等研究的基础上考虑了原位碳释放过程和硝化作用,揭示了大坝释水时河水位波动条件下潜流带氮的迁移规律,但却忽略了河床和含水层非均质性在潜流带氮循环中起到的影响,且研究中假设含水层是均质的。当前已有大量的野外观测和研究表明,河床和含水层的非均质性会影响潜流交换[18-19]、溶质停留时间[20-22]和氮的去除程度[23-25]。一般来说,在非均质场中渗透系数较大的区域,溶质滞留时间较短且氧气浓度较高[26],不利于反硝化作用的发生;渗透系数较小的区域,溶质停留时间较长,有利于消耗氧气形成厌氧环境,利于硝酸盐的去除[27-29]。Sawyer等[30]采用数值模拟手段研究发现,河床沉积物非均质性影响进入潜流带内的水流运动轨迹和在潜流带内的停留时间;Bardini等[31]采用数值模拟方法在研究河床形态对氮素迁移转化的影响时发现,泥沙渗透性和河流流速直接影响了养分的供应和在河床中的停留时间,控制了沙丘下的反应速率;Pescimoro等[32]在研究地梯度蜿蜒河流中具有砂和淤泥二元分布沉积物中的潜流交换时发现,虽然更高淤泥占比的河床表现出较低的潜流交换率,但是沿流动路径上硝酸盐的去除效率更高;张强[33]建立了半切割河床式潜流带模型,并将均质潜流带模型和低非均质性潜流带模型进行了对比,发现潜流带非均质性降低了潜流交换量和硝酸盐削减率;Wallace等[34]采用数值模拟手段研究了沉积物非均质性对潮汐河流附近氮动态的影响,发现随着含砂量的增加,河水中去除的硝酸盐数量增加。

尽管目前关于河水位波动和渗透系数非均质对潜流带生物地球化学过程影响的研究已有很多,但大多研究往往忽略了河水位波动过程中氧气溶解和原位碳释放溶解有机碳的过程,这可能会导致非均质潜流带模型低估了潜流带中溶质供给和反硝化作用。为了进一步探究潜流带氮循环过程,本研究建立潜流带二维流动和反应传输模型,探究河水位波动下渗透系数非均质性与潜流带内氮素迁移转化过程的空间关系,了解渗透系数非均质性对潜流交换的影响、溶质在空间上迁移转化的差异和硝酸盐去除率的变化。

1 非均质潜流带模型的建立

1.1 概念模型建立

使用COMSOL Multiphysics构建二维潜流带概念模型,用于模拟河流波动引起的潜流带内氮素迁移转化过程。模型几何形状取自Hornsby Bend站点上用全站仪测量的从河漫滩到河道中心的河岸高程轮廓,具体见图1。

注:透水层外部高度为河水位高度;开放边界外部浓度为河水浓度和以气泡形式捕获的空气浓度。图1 二维潜流带氮循环过程示意图(改自Shuai等[17])Fig.1 Schematic diagram of the nitrogen cycle process in a two-dimensional hyporheic zone (adapted from Shuai et al.[17])

1.2 数学模型建立

本研究采用Richard方程和对流-弥散-反应方程刻画潜流带内水流和溶质的迁移过程,其控制方程如下:

(1)

(2)

式中:p为压力(Pa);g为重力加速度(m/s2);ρ为流体密度(kg/m3);μ为流体动力黏度(Pa·s);Cm为比水容量(1/m);Se为有效饱和度(无量纲);S为储水系数(1/Pa);ks为饱和水力渗透率(m2);kr为相对渗透率(无量纲),是有效饱和度的函数;t为时间(s);z为位置高度(m);Qm为源汇项[kg/(m3·s)];θ为含水率(m3/m3);Ci为物质i的浓度(mg/L);u为达西流速(m/s);D为水动力弥散系数(m2/s);Ri为物质i的反应速率[mg/(L·d)]。

有氧呼吸作用 CH2O+O2→CO2+H2O

原位有机物或总有机碳的溶解可能导致DOC的输入,该过程通过一级传质过程来模拟[17,37-40]:

(3)

式中:CDOC和CTOC分别为沉积物中DOC和总有机碳(TOC)的浓度(mg/L);α为一阶传质系数(1/h);Kd为DOC分布系数(L/kg)。

有氧呼吸作用、硝化作用和反硝化作用由多重Monod动力学描述[41]。存在DO时,使用非竞争性来抑制反硝化作用[42]。此外,系统中的氧通常用于能量优先级更高的反应,如有氧呼吸;而剩余的氧则作为具有较低能量优先级的二级溶质被消耗,如反硝化反应。而溶解氧需求的分布系数yO2是基于氧化和硝化反应的吉布斯自由能计算的。综上所述,各化学反应项表示如下:

(4)

(5)

(6)

(7)

1.3 模型网格剖分与边界设置

基于概念模型和数学模型,首先将潜流带模型离散化。为了提高模型模拟精度和减少计算量,采用三角形网格进行剖分(图2),上边界附近最小空间步长为0.1 m,远离上边界最大空间步长为0.3 m,将模型离散为10 184个网格。模型模拟总时长为24 h,时间步长为1 h。

对于水流模型,河道中心剖面与底部边界处理为隔水边界;右侧边界和上边界分别设置为定水头边界和透水层边界,不考虑降雨或蒸发情况(具体设置参考Cardenas等[43])。河水位波动近似表示成振幅为0.5 m、周期为24 h的正弦函数。

对于溶质运移模型,河道中心剖面与底部边界处理为零通量边界;右侧边界设置为流出边界;上边界设置为开放边界,在这里允许发生溶质的流入和流出,且外部浓度等于河流中溶质浓度。假定在地下水水位上升时,与大气相连的非饱和带中的DO处于饱和状态。

表1 水流场和溶质运移场主要参数设置

2 结果与讨论

2.1 潜流交换通量和流场分布情况

渗透系数的非均质性影响河流-潜流带的水动力交换和流场分布,一个24 h周期内河水位波动和潜流交换通量,见图3。

注:正值表示含水层向河流补给(出渗),负值表示河水向含水层补给(入渗);平均河水位和零通量则分别用红色和蓝色虚线表示。图3 一个24 h周期内河水位波动(蓝色实线)和潜流 交换通量(点线)Fig.3 Fluctuations in river level (solid blue line) and hyporheic exchange flux (dotted line) in a 24-hour period

由图3可知:河水位上升时,河水向含水层入渗(负值),河水位下降至平均河水位以下后,水由河水入渗含水层(负值)转变为流出含水层(正值)向河水补给,并在18 h时河水位达到最低;18 h后,河水位虽然有所上升,但此时河水位仍低于平均河水位,含水层入渗补给河水;河流-含水层交界面上的入渗和出渗通量随着河水位波动而上升和下降,在河段河水位达到最大值前入渗流量达到峰值,而出渗通量在河段河水位达到最小值前达到峰值,这一趋势与 Musial等[9]受潮汐影响的近溪流实地观察结果和Shuai等[17]的模拟结果一致。

注:黑色箭头为流线。图4 不同方差下潜流带渗透系数场和流场分布图Fig.4 Diagram of permeability coefficient fields and flow fields of hyporheic zones under different variance values

由图4可知,随着潜流带含水层非均质性的增强,流线越曲折,此外,含水层的非均质性还影响潜流带的交换速率,含水层的非均质性越强,潜流带交界面处的等水位线向两侧延申越远,潜流带交换范围越广,反应越剧烈。

2.2 溶质浓度和生物地球化学反应速率

图5 不同方差下潜流带中溶质运移过程图Fig.5 Solute transport process in hyporheic zones with different variance values

图6 不同方差下潜流带中观测点处溶质浓度和化学反应速率对比图Fig.6 Comparison between solute concentration and chemical reaction rate at observation points in hyporheic zone under different variance values

硝化作用受溶解氧控制,主要发生在河流-含水层界面附近的富氧区域[图5(d)];反硝化作用发生在厌氧区域,主要沿着有氧呼吸区域外围以带状的锋面形式分布[图5(e)],越远离高氧区,反硝化反应速率就越高。渗透系数的非均质性使潜流带内高渗透层增多,溶质的渗透深度增加,氨氮、硝酸盐等物质向潜流带深处运移,这在一定程度上了改变了硝化和反硝化带的形态,扩大了反应的作用区。与潜流带内溶质浓度的空间分布形状相类似,随着含水层非均质性的增强,潜流带内高渗透层增加,硝化作用和反硝化作用区分布形状越不规则,边缘越“曲折”。

2.3 硝酸盐去除率

为了定量分析潜流带对硝酸盐的去除情况,Shuai等[17]假定硝酸盐去除效率是含水层中去除的硝酸盐总量与进入含水层中硝酸盐总量的比值NRE,其可表示为

其中,反硝化作用去除的含水层中硝酸盐量MDN和经河流-含水层界面进入含水层中的硝酸盐总量Min分别表示如下:

渗透系数的非均质性影响了潜流带内硝酸盐的去除,不同方差下潜流带内硝酸盐的去除情况见表2。

表2 不同方差下潜流带内硝酸盐去除情况

由表2可知:随着含水层非均质性程度的增加,硝化、反硝化反应速率和溶质渗透深度均有所增加。总得来说,一个周期内河流向潜流带输入的硝酸盐总量和经反硝化作用去除的硝酸盐总量也随之增加,潜流带内硝酸盐去除率随着渗透系数非均质性程度的增加而增大。

3 结 论

本研究利用数值模拟手段研究了渗透系数非均质性对潜流带氮素迁移转化和硝酸盐去除的影响。主要得到结论如下:

1) 潜流内渗透系数非均质性程度的增加加快了潜流交换并提供了更高的反应速率。随着非均质性程度的增加,潜流带渗流面上入渗率和出渗率也越大,河水-含水层间的水力交换越剧烈。

2) 潜流内渗透系数非均质性影响了溶质在含水层中的运移情况,控制了潜流带内的反应动力学,改变了溶质分布形状和硝化、反硝化作用的发生区域,随着非均质性的加强,溶质浓度和反应区在空间的分布越不规则,分布区边缘越“曲折”。

3) 渗透系数的非均质性影响潜流带去除硝酸盐的能力,更高的渗透系数非均质性表现出更好的硝酸盐去除效果,硝酸盐去除率也有所增加。