利用HYDRUS-2D软件模拟污染事故后三氮污染物的迁移转化规律*

尹芝华 杜青青 翟远征 杨 洁 郭永丽 左 锐#

(1.北京师范大学水科学研究院，地下水污染控制与修复教育部工程研究中心，北京 100875;2.中国地质科学院岩溶地质研究所，国土资源部/广西岩溶动力学重点实验室，广西 桂林 541004)

利用HYDRUS-2D软件模拟污染事故后三氮污染物的迁移转化规律*

尹芝华1杜青青1翟远征1杨 洁1郭永丽2左 锐1#

(1.北京师范大学水科学研究院，地下水污染控制与修复教育部工程研究中心，北京 100875;2.中国地质科学院岩溶地质研究所，国土资源部/广西岩溶动力学重点实验室，广西 桂林 541004)

以某污水处理池泄露污染水环境为例，利用HYDRUS-2D软件构建土壤水分运动和溶质运移模型，模拟三氮在该场地非饱和带垂向以及向下游地表水体的迁移转化过程，为定量预测识别该污水泄漏场地对下伏含水层和下游河流的影响提供技术参考。结果表明：(1)三氮在非饱和带垂向100cm以内明显富集，但富集规律各异；(2)氨氮和亚硝态氮在整个模拟期内并未迁移至潜水面，对下游河流的影响微乎其微，最终污染下伏含水层和下游河流的是硝态氮；(3)硝态氮在污水泄漏事件发生后的第1 095天以0.024 0mg/cm3(即24mg/L)到达下伏潜水面，超过《地下水质量标准》(GB/T14848—93)中Ⅲ类限值(20mg/L)；(4)硝态氮在污水泄漏事件发生后的第938天以0.468 5mg/cm3(即468.5mg/L)的峰值到达下游河流，超过《地表水环境质量标准》(GB3838—2002)中规定的Ⅴ类限值(2.0mg/L)。

三氮 非饱和带 迁移转化 数值模拟HYDRUS-2D

近年来，污染场地作为一个世界性环境问题，对人类健康和生态环境产生了严重威胁，并显著影响了土地和地下水资源的安全使用[1]。污染场地是由多种介质和种类繁多的污染物组成的复合系统，其中三氮在非饱和带中受复杂的物理化学生物过程影响，是污染场地地下水中典型的无机污染物。国土资源部发布的《2015中国国土资源公报》显示，全国202个地市级行政区5 118个地下水水质监测点中，42.5%的点位水质呈较差或极差状态，主要超标组分为三氮、总硬度、溶解性总固体、铁、锰等。位于地下水位以上的非饱和带是控制地下水污染的关键区域，非饱和带作为三氮由地表向地下传输的过渡带，对迁移其中的三氮进行滞留阻截和消散，是保护地下水的天然屏障[2]。因此，研究三氮在非饱和带中的迁移转化过程是控制地下水污染的基础，对保护水环境具有重要意义。现有研究主要针对三氮通过非饱和带对地下水[3-6]和地表水体[7-11]的污染过程及程度等方面进行探讨，对三氮在非饱和带迁移转化后对地下水和地表水体影响的联合研究较少。

对于非饱和带中的水分运动和溶质运移过程研究，由于空间和时间上的强烈变异性，需要采用大量的测量工作且较为复杂[12-14]。随着计算机技术的迅速发展，采用数值模拟方法模拟非饱和带土壤水分运动和溶质运移越来越受人们关注，且是目前解决复杂条件下土壤水分运动问题的有效方法[15]。本研究以某污水处理池泄露污染水环境为例，利用HYDRUS-2D软件对非饱和带构建水流运动和溶质运移模型，模拟污水处理池渗漏液中的特征污染物三氮在非饱和带垂向以及向下游地表水体的迁移转化过程，并预测三氮到达下伏潜水含水层和下游河流的时间和稳定浓度值，为定量预测识别该污水泄漏场地对下伏含水层和下游河流的影响提供依据。

1 研究区概况

污水处理池位于某沙漠边缘，长8 m，宽6 m，深2 m，污水成分主要为氨氮，由于年久失修可能会发生泄漏情况。本次假定未来某一时刻泄漏事件发生，且3年后被发现，模拟污水成分对下伏含水层及下游河流的影响程度，单位面积渗漏量根据已有资料确定为0.003 6 m3/d。

污水处理池所在场地整体地势西高东低，场区西侧为山区，东部约30 m处有河流自北向南流过，该河流对于附近居民农业方面有重要意义。场地所在区域属温带大陆性干旱气候区，日照充足，昼夜温差大，年降水量和蒸发量很少。经过实地调查，场地分布有灰漠土，边缘有沙化灰漠土，绝大多数地表为农田。主要的土种有立土、平土和漏沙土，类型以沙壤土为主，渗透系数18.76 cm/d。地下水类型为第四系松散层孔隙水。由于当地荒无人烟，地下水流场处于天然状态，水位埋深约5 m。

2 模型建立

2.1 场地条件概化

根据场地水文地质勘查结果，下伏含水层接近非饱和带底部。污水处理池西侧为山区，东侧为河流。存在河流补给地下水的微弱水力联系。水流模型可概化为非均质各向同性多孔介质非饱和剖面二维非稳定流，上边界为地表，下边界为潜水面，污水处理池泄漏点距潜水面约380 mm，研究区左右两侧为定压力水头。土壤溶质运移模型根据多孔介质溶质运移理论，考虑三氮之间的转化。

经过工程勘察后，选取典型土壤样品在实验室测定土壤特征参数，以确定后期模拟相关参数。本次模拟选择渗漏液中浓度及影响较大的三氮作为特征污染物，其中氨氮初始泄漏水平根据污水处理池污水成分检测结果所定，为1.152 0 mg/cm3。

2.2 数值模型

目前，非饱和带多孔介质中的水分运动和溶质运移过程通常由Richards方程和对流-弥散方程(CDE)来刻画[16]。运用Richards方程描述水分运动时，通常不考虑液体流动过程中气体和热梯度的影响。Richards方程中土壤水力特性参数由目前应用最广泛的van Genuchten-Mualem模型[17]计算，且不考虑水流运动的滞后现象。

运用CDE[18]描述三氮迁移转化过程时，假定三氮之间的转化过程是一个全程硝化反硝化过程，两个反应独立进行，即氨氮通过亚硝化作用生成亚硝态氮，亚硝态氮通过硝化反应生成硝态氮，反硝化过程中对硝态氮的代谢全部为异化代谢，即全部被还原成气态氮，而不发生同化代谢生成氨氮和亚硝态氮。

2.3 初始条件和边界条件

根据研究区域周边实际情况，选择定压力水头边界，模型左、右侧边界水头分别设置为400、180 cm，保持线性平衡分布。对于水分运动上边界，河流周围定为定压力水头边界，设置为80 cm，河流左侧设置为渗出面边界，模型下边界概化为自由排水边界。

溶质运移侧向边界与水分运动侧向边界相一致。垂向边界设置为第三类边界，在河流左侧设置浓度边界，边界浓度根据野外钻孔中所取土样的浸出试验结果取平均值给定，氨氮质量浓度为1.152 0 mg/cm3，下边界设置为浓度零梯度边界。

由于氮的转化过程十分复杂，建模时假设三氮的转化过程遵循氨氮转化为亚硝态氮再转化为硝态氮，忽略土壤中有机氮的矿化作用和土壤颗粒对亚硝态氮和硝态氮的吸附。同时，假定氨氮和亚硝态氮的硝化作用和硝态氮的反硝化作用均符合一级动力学反应规律。本次模拟假定非饱和带中氨氮、亚硝态氮、硝酸盐氮的初始质量浓度均为零，即假定非饱和带在污水处理池泄漏前未接受氮源补给。

表1 土壤水力参数

表2 溶质运移模型参数

注：1)指在纯水中的分子扩散系数。

2.4 参数率定

(1) 水分运动参数率定方面，土壤水力参数包括土壤的残余含水率(θr，cm3/cm3)、土壤的饱和含水率(θs，cm3/cm3)、进气值倒数(α，m-1)、孔径分布指数(n)、渗透系数(Ks，cm/d)、土壤介质孔隙连通性能参数(l)，其中Ks为实测值，其他土壤水力参数根据HYDRUS-2D软件内置Rosetta模型，通过输入实测土壤粒径比例预测参数初值，经过调整后的van Genuchten-Mualem公式中各个土壤水力参数值见表1。

(2) 溶质运移参数率定方面，三氮的吸附主要表现为线性吸附，吸附能力高低取决于分配系数的大小，一般来说，分配系数越大越易吸附而越不易迁移。沙壤土的吸附分配系数通过以下公式计算：

Kd=1.34+0.083CL

(1)

式中：Kd为吸附分配系数，cm3/mg；CL为黏粒质量分数，%，数值从《第四纪松散沉积物分类命名表》[19]中查得。

弥散度具有很大的尺度效应，实验室与野外现场所测数值可以相差数个数量级[20]。王全九等[21]研究发现，可以根据非饱和带土壤饱和导水率推求弥散度。其中，横向弥散度大约为纵向弥散度的1/10～1/3[22]。土壤容重参考不同岩性土壤密度经验值[23]确定。主要溶质运移参数如一级反应速率常数等在以上数据基础上确定，具体方法为：在软件中反演程序中输入不同时间和土壤深度的溶质浓度实测值，然后通过反演对参数进行优化，从而得到参数值，率定后的参数见表2。

2.5 空间离散和时间离散

数值求解过程中，采用Galerkin有限元法对土壤剖面进行空间离散，采用隐式差分格式进行时间离散。根据研究区的边界条件和水动力条件，在整个模拟区域内，把模型剖分成1 513个节点和2 859个三角形单元网格。为更加直观地表现三氮垂向变化，对纵向刻度做扩大处理。对于时间步长，选择模型内部推荐步长，即初始时间步长为110-4d，最小时间步长为110-5d。模拟总时长为5年(1 825 d)，其中前3年(1 095 d)为污水泄漏期。整个模拟区共选择1个观测孔M1，其在模拟区的位置和模型网格剖分见图1。

图1 模型网格剖分Fig.1 Grid splitting map of model

3 结果与讨论

3.1 模型验证

由于场地现场验证困难大，为检验HYDRUS-2D软件模拟效果的可靠性，在研究区采集土壤和污水开展室内土柱试验，并使用HYDRUS-2D软件对土柱内土壤中氨氮的迁移转化过程进行模拟，用试验结果对模拟结果进行验证。土柱试验设计模拟剖面深度为65 cm，将24 h后的模拟值和土柱试验24 h后的实测值进行对比，结果见图2。由图2可见，拟合值与实测值拟合效果较好，说明所建数值模型比较可靠，可用于模拟预测。

3.2 三氮在非饱和带中的演化规律预测

在上述建模与参数识别、模型验证的基础上，经过数值模拟过程得到三氮在非饱和带中的迁移转化规律，其中模拟剖面上氨氮、亚硝态氮、硝态氮的污染晕分别见图3至图5。根据模拟结果，在非饱和带中，氨氮、亚硝态氮在模拟期结束前已停止运移且质量浓度减为0 mg/cm3，硝态氮在整个模拟期内都存在。

图2 氨氮质量浓度实测值与模拟值的对比Fig.2 Comparison of measured and simulated values of ammonia nitrogen

3.2.1 三氮在非饱和带中的垂向迁移转化

在污染点垂直方向上，三氮在第50、365、730、1 095天的分布情况见图6。总体看来，沙壤土介质中三氮在非饱和带垂向100 cm以内明显富集，但迁移规律略有差异。氨氮从地表进入土层后，随水分向下迁移过程中浓度迅速降低，但迁移速度较慢。整个模拟期氨氮的迁移主要发生在地表以下300 cm以内。在第50天，氨氮在埋深0～180 cm的剖面上从1.152 mg/cm3迅速降到0 mg/cm3。第365天氨氮在土壤表层开始出现累积，但累积量并不大，大约在地下10 cm处累积量达到1.171 0 mg/cm3，迁移深度达到地下280 cm深处。第730天最大累积量出现在地下约16 cm处，达到1.197 0 mg/cm3，迁移深度达到地下340 cm深处。第1 095天氨氮在表层逐步累积，向下迁移现象不明显，最大累积量出现在地下32 cm深处。亚硝态氮随水分向下运移速度较快，但到第730天后较少向下迁移。整个模拟期亚硝态氮的迁移主要发生在地表以下375 cm以内。在第50天，亚硝态氮在土壤表面出现了累积现象，最大累积量出现在约19 cm深处，约为1.042 0 mg/cm3，并向下迁移了212 cm。第365天亚硝态氮在地下约26 cm深处达到最大累积量1.051 0 mg/cm3，迁移深度达到地下330 cm深处。第730天最大累积量出现在地下约36 cm深处，达1.061 0 mg/cm3，迁移深度达到地下352 cm深处。第1 095天亚硝态氮最大累积量出现在地下约42 cm处，达1.064 0 mg/cm3，迁移深度达到地下375 cm深处，但并未迁移至潜水面。硝态氮随水分向下运移速度较快，且在土壤表层同样有累积现象。在第50天，硝态氮最大累积量出现在约23 cm深处，为1.036 0 mg/cm3，并向下迁移190 cm。第365天硝态氮在地下约38 cm处达到最大累积量1.063 0 mg/cm3，迁移深度达到地下280 cm深处。第730天最大累积量出现在地下约42 cm深处，达到1.074 0 mg/cm3，迁移深度达到地下326 cm深处。第1 095天硝态氮最大累积量出现在地下约44 cm深处，达到1.080 0 mg/cm3，且随水流迁移至潜水面，最终进入潜水面的硝态氮为0.024 0 mg/cm3，即24 mg/L，超过《地下水质量标准》(GB/T 14848—93)中硝态氮的Ⅲ类限值(20 mg/L)；第1 095天后，由于泄露事件停止，氨氮和亚硝态氮含量急剧下降，硝化反应作用下降，三氮累积量逐渐降低，小于峰值浓度。

图3 模拟剖面上氨氮的污染晕变化Fig.3 The variation of ammonia nitrogen pollution halo in the simulated profile

图4 模拟剖面上亚硝态氮的污染晕变化Fig.4 The variation of nitrite pollution halo in the simulated profile

图5 模拟剖面上硝态氮的污染晕变化Fig.5 The variation of nitrate pollution halo in the simulated profile

图6 不同时间的三氮垂向分布Fig.6 Vertical distribution of three nitrogens concentration in different time

总体来看，氨氮进入土壤后，随深度增加，其迁移能力减弱，考虑到土壤对氨氮的吸附及氨氮自身的转化作用等因素，氨氮进入土层后浓度迅速减少，同时在土壤表层出现累积现象，说明土壤表层对氨氮的吸附已趋于饱和。相似地，亚硝态氮和硝态氮也在土壤表层出现累积现象，但土壤颗粒对亚硝态氮和硝态氮几乎不吸附，而且沙壤土岩性颗粒粗，透水性好，使得亚硝态氮和硝态氮在垂向上迁移速度较快。随着水样的不断垂向渗透，非饱和带毛细效果减弱，土壤孔隙中自由水含量减小，饱和度降低，氧化性增强，促进硝化反应的进行，加上亚硝态氮本身的不稳定性，极易通过硝化反应转为硝态氮，故最后通过非饱和带进入地下水的三氮组分主要为硝态氮。同时，在迁移过程中，土层中的含氧量随深度逐渐降低，硝化作用由强变弱，到达一定非饱和带深度后硝化作用几乎为零。所以垂向上，氨氮、亚硝态氮和硝态氮浓度分布上总体上高下低，在一定深度以下浓度降至零。

3.2.2 三氮向河流的迁移转化

M1观测孔位于地表水和地下水交错带的核心地带，其三氮质量浓度的时间变化见图7。由图7可见，氨氮在非饱和带经过迁移转化作用后，进入河流的浓度微乎其微，可忽略不计，亚硝态氮迁移趋势与氨氮相似。硝态氮浓度整体上呈现出先迅速升高至稳定并维持一段时间后再迅速降低的发展规律。硝态氮从第490天左右进入河流，在第938天达到峰值，为0.468 5 mg/cm3(即468.5 mg/L)，超过《地表水环境质量标准》(GB 3838—2002)中规定的TN Ⅴ类标准限值(2.0 mg/L)。第1 600天后河流中的硝态氮质量浓度急剧降低，在第1 825天时，进入河流的硝态氮的质量浓度为0.058 0 mg/cm3。可见，氨氮和亚硝态氮在非饱和带中的迁移与转化过程非常迅速。同时，由于沙壤土土质较为疏松，岩性颗粒粗，空隙稍大，其中含氧量较高，促进硝化反应的进行，故硝态氮含量不断增高，直至稳定状态，泄漏停止后，氮源减少，硝态氮在反硝化作用下，浓度逐渐降低。

图7 M1观测孔中三氮质量浓度的时间变化Fig.7 Time variation of three nitrogens concentration in M1 observation well

4 结 论

(1) 模拟场地的非饱和带介质以壤土为主，厚度约500 cm，是氮污染负荷的有效缓冲区，对下伏含水层和下游河流能起到一定的保护作用，但是由于非饱和带对亚硝态氮和硝态氮的吸附能力相对较弱，导致两者的下渗量相对较大。而亚硝态氮不稳定，大部分经过硝化过程转化为相对稳定的硝态氮，因此硝态氮为主要污染物。

(2) 模拟期内三氮在非饱和带垂向100 cm以内明显富集，氨氮和亚硝态氮并未到达潜水面，硝态氮在第1 095天以0.024 0 mg/cm3(即24 mg/L)到达潜水面，超过GB/T 14848—93中规定的Ⅲ类限值(20 mg/L)。

(3) 模拟期内氨氮和亚硝态氮进入河流的浓度微乎其微，可以忽略不计。硝态氮在第490天到达地表水体，并在第938天达到稳定质量浓度0.468 5 mg/cm3(即468.5 mg/L)，远超GB 3838—2002中规定的Ⅴ类限值(2.0 mg/L)。

[1] 廖晓勇,陶欢,阎秀兰,等.污染场地修复决策支持系统的几个关键问题探讨[J].环境科学,2014,35(4):1576-1585.

[2] 郭敏丽,王金生,杨志兵,等.非正规垃圾填埋场包气带介质的污染物阻滞能力研究[J].环境污染与防治,2012,34(1):24-26.

[3] 陈英旭,梁新强.氮磷在农田土壤中的迁移转化规律及其对水环境质量的影响[M].北京:科学出版社,2012.

[4] 黄亮,吴乃成,唐涛,等.水生植物对富营养化水系统中氮、磷的富集与转移[J].中国环境科学,2010,30(1):1-6.

[5] 刘海军,李艳,张睿昊,等.北京市集约化种植土壤硝态氮分布和迁移速率研究[J].北京师范大学学报(自然科学版),2013,49(2):266-270.

[6] GU B,GE Y,CHANG S X,et al.Nitrate in groundwater of China:sources and driving forces[J].Global Environmental Change,2013,23(5):1112-1121.

[7] 左锐,马丁山,刘丽,等.傍河水源地氨氮分布特征及其影响因素分析[J].北京师范大学学报(自然科学版),2015,51(6):626-630.

[8] 赵丽,张韵,张丹,等.重庆西部红层浅层地下水中“三氮”污染现状及影响因素分析[J].环境污染与防治,2016,38(6):46-53.

[9] 陈海洋,滕彦国,王金生,等.晋江流域非点源氮磷负荷及污染源解析[J].农业工程学报,2012,28(5):213-219.

[10] 李远,常学礼,孙朋,等.降水波动背景下嫩江干流氨氮污染变化趋势分析[J].环境污染与防治,2013,35(5):1-5.

[11] 傅国伟,程声通.水污染控制规划[M].北京:清华大学出版社,1985.

[12] 郭劲松,李胜海,龙腾锐.水质模型及其应用研究进展[J].土木建筑与环境工程,2002,24(2):109-115.

[13] 丁素玲.非均质包气带中“三氮”迁移转化数值模拟[D].北京：中国地质大学(北京),2012.

[14] 董晓华,刘潇钧,彭涛,等.利用Hydrus-1D模拟分层土壤剖面的水流运动[J].安徽农业科学,2016,44(4):29-31.

[15] 邵明安,王全九,黄明斌.土壤物理学[M].北京:高等教育出版社,2006.

[16] 杨金中,蔡树英,王旭升.地下水运动数学模型[M].北京:科学出版社,2009.

[17] VAN GENUCHTEN M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Science Society of America Journal,1980,44(5):892-898.

[18] 王小丹,凤蔚,王文科,等.基于HYDRUS-1D模型模拟关中盆地氮在包气带中的迁移转化规律[J].地质调查与研究,2015,38(4):291-298.

[19] 何江涛，马文洁，张听.基于过程单元模型参数替代防污性能评价方法研究——以北京市平原区为例[M].北京：地质出版社，2012.

[20] 许秀元,陈同斌.土壤中溶质运移模拟的理论与应用[J].地理研究,1998,17(1):99-106.

[21] 王全九,邵明安,郑纪勇.土壤中水分运动与溶质迁移[M].北京:中国水利水电出版社,2007.

[22] 郑春苗.地下水污染物迁移模拟[M].2版.北京:高等教育出版社,2009.

[23] 唐大雄,刘佑荣,张文殊,等.工程岩土学[M].2版.北京:地质出版社,1999.

SimulationandpredictionofleakageincidentonmigrationandtransformationofthreenitrogensbasedonHYDRUS-2Dsoftware

YINZhihua1,DUQingqing1,ZHAIYuanzheng1,YANGJie1,GUOYongli2,ZUORui1.

(1.CollegeofWaterSciences,BeijingNormalUniversity,EngineeringResearchCenterofGroundwaterPollutionControlandRemediation,MinistryofEducation,Beijing100875;2.InstituteofKarstGeology,ChineseAcadamyofGeologicalScience，KarstDynamicsLaboratory,MinistryofLandandResources&Guangxi,GuilinGuangxi541004)

A sewage treatment tank leakage of water pollution was taken as an example to construct the soil water migration and solute transport model by using HYDRUS-2D software. The migration process of the three nitrogens in the unsaturated zone and the downstream surface water were simulated,providing the technical reference of quantitative prediction to identify the impact of the wastewater spill site on the underlying groundwater and downstream rivers. The results show that: (1) three nitrogens were enriching within 100 cm of the unsaturated zone,but the enrichment regulations were different. (2) Ammonia nitrogen and nitrite nitrogen did not migrate to the diving surface during the whole simulation period and their impacts on the downstream river were negligible. The final contamination entering the lower aquifer and downstream river was the nitrate. (3) Nitrate nitrogen reached the lower diving surface at a stable peak concentration of 0.024 0 mg/cm3(24 mg/L) on the 1 095th day after the occurrence of the sewage spill,exceeding the Ⅲ limit of nitrate nitrogen (20 mg/L) in the “Groundwater quality standard” (GB/T 14848-93). (4) Nitrate nitrogen reached the downstream river at a peak concentration of 0.468 5 mg/cm3(468.5 mg/L) on the 938th day after the occurrence of the sewage spill and simultaneously the total nitrogen content exceeds the total nitrogen Ⅴ limit (2.0 mg/L) stipulated in the “Surface water environmental quality standard” (GB 3838-2002).

three nitrogens; unsaturated zone; migration and transformation; numerical simulation; HYDRUS-2D

10.15985/j.cnki.1001-3865.2017.10.006

尹芝华，男，1991年生，硕士研究生，研究方向为非饱和-饱和带污染迁移的数值模拟。#

。

*国家水体污染控制与治理科技重大专项(No.2014ZX07201-010)；岩溶动力学重点实验室基金资助项目(No.KDL201505)。

2017-03-06)