牧区河岸潜流带硝酸盐氮和氨氮浓度对水文过程的响应机制

薛栋元， 胡海珠，2， 张锦宁， 任嘉伟

(1.内蒙古大学生态与环境学院，内蒙古河流与湖泊生态重点实验室，内蒙古 呼和浩特 010020；2.蒙古高原生态学与资源利用教育部重点实验室，内蒙古 呼和浩特 010020)

氮素引起的富营养化是全球范围关注的水环境问题［1］。农业生态系统中合成肥料的过量施用和畜牧养殖业中畜禽粪尿的大量排出是氮素面源污染的重要来源。大气、地表和土壤中的过量氮素在水文过程的驱动下进入水体，从而造成水环境污染，硝酸盐（）和氨盐（）是水体中最主要的氮素形态［2］。河岸潜流带作为地表水和地下水的混合区域，其中的水量交换促进了物质和能量的交换，并形成了明显的氧化还原梯度和多样的生物群落，因此被普遍认为是营养物质转化，特别是氮素循环的热点区域［3-6］。厘清两种氮素在河岸潜流带的迁移转化和来源是氮素污染负荷控制的关键。

有研究表明，夏季汛期是氮素转化的热点时期，降雨径流事件是流域氮素面源污染的主要驱动因子，刺激了氮循环的周转［7］。相关研究通过分析污染物浓度与降雨、流量的动态关系，利用滞后现象（hysteresis）解释氮素的来源和途径，即氮素是来自于涨水期通过地表径流快速进入河流或是在落水期由地下水排泄至河流［8-9］，但是仅依靠河水的流量和氮素浓度监测及滞后分析不足以准确识别氮素的迁移转化，需要结合氮素及相关环境因子在河岸带的分布进一步明晰其周转过程。

大多数研究表明，河岸潜流带是硝酸盐的汇，其中发生的硝化、反硝化、厌氧氨氧化等过程取决于潜流带中的氮素供应、水化学特征以及水文条件（例如洪水特征、地下水位）等［10-11］。目前，研究者对灌溉农业区附近的河流潜流带开展了广泛研究，例如，Darwiche-Criado等［12-14］发现，在丰水期与枯水期的来源不同，且对水文过程的响应不同。潘俊等［15-17］发现了潜流带的氧化还原条件影响了氮素的生物地球转化过程，反硝化、厌氧氨氧化转化强度在一定范围内存在空间差异。但是，目前国内研究对牧区河岸带的氮素迁移转化规律认识仍有不足，大多研究仅关注在反硝化作用下的迁移转化模拟［18-19］。基于此，本研究基于汛期高频率水位、水体氮素浓度及其相关环境因子的监测，采用模拟软件FEFLOW［20］构建河岸潜流带水流-氮素迁移转化模型。旨在分析河流及河岸潜流带中氮素动态及其对水文过程的响应，探究河岸潜流带中两种主要氮素溶质的来源和运移机制，以期为河岸带管理及流域氮素污染治理提供借鉴与帮助。

1 材料与方法

1.1 研究区概况

研究区域位于内蒙古高原东部的锡林郭勒草原，地处我国草原从东部半湿润草甸草原区向西北半干旱典型草原区过渡地带（图1）。锡林河是锡林郭勒草原的母亲河，全长268 km，从东南向西北流动。锡林河流域面积10542 km2，年降水量为250～400 mm，全年最大降水量出现在夏季汛期的7—8月，多年平均蒸发量1900 mm［21］。地下水以松散岩孔隙水为主，富水性较强，上层潜水含水层以上更新统的细砂为主［22］。锡林郭勒草原是我国北方重要的牧场。

图1 研究区及研究河段监测井布设Fig.1 Overview of the study area and layout of the monitoring wells in the studied river section

1.2 原位定点监测试验

于2020年夏季汛期，在位于中国科学院内蒙古草原生态系统定位研究站附近的锡林河上游河段开展野外定点监测试验。监测时段为6 月15 日、7月2 日至8 月14 日、9 月20 日。研究河段总长110 m，河宽平均8.5 m。监测期内降雨量共180.3 mm，其中两次较大降雨，分别为7 月11—13 日连续3 d降雨共47 mm 和8 月8 日的降雨64 mm。监测期内河水深0.3～0.7 m，河岸地下水埋深为0.5～1.1 m。分别在河道及河岸共布设9个监测井（图1），河道内间隔5 m 均匀布设3 个监测井(C1～C3)，河道左岸（1L～3L）和右岸（1R～3R）分别布设3个监测井，位于距离河道7～10 m 位置。河岸监测管深度为地表以下2 m，河道监测管深度为河床表面以下1 m。选用直径为50 mm的PVC管作为地下水监测井管，在距离管底10～15 cm 的位置均匀钻孔，并用400 目尼龙网包裹，防止淤泥堵塞。

监测前期，采用全站仪测量研究区地形。监测期间，采用水位计及水尺逐日监测地下水位及河水位。使用流速仪测量河水流速，根据河宽换算流量。每隔5～7 d 采用便携式水质分析仪（Multi 3630 IDS SETG）监测河水和地下水的温度、溶解氧（Dissolved Oxygen，DO）、pH 等水质参数。同时采集河水和地下水水样，置入100 mL 棕色采样瓶，阴凉且避光保存，防止发生生物化学转化。将采集水样在实验室通过直径0.4 μm滤膜过滤，根据生态环境部推荐的测定方法，使用Vario TOC分析仪测定溶解有机碳（Dissolved Organic Carbon，DOC）浓度［23］。使用紫外分光光度计（SHIMADZU UV-2600）测定无机氮浓度：采用纳氏试剂分光光度法测定［24］；采用紫外分光光度法测定［25］。

1.3 河岸潜流带水流和溶质运移模型构建

研究区河岸潜流带通过降水和河流渗透得到补给，河水和地下水的交换取决于河流与地下水位之间的水力梯度和沉积物的渗透性能。放牧牛群排泄的粪尿是研究区主要的氮素来源，在水文和生物地球化学的共同作用下可能发生硝化、反硝化等反应过程。本研究基于FEFLOW 建立了研究河段河岸潜流带的三维水流和氮素溶质反应运移数值模型（图2），根据潜流带内水流过程和氮素迁移转化过程，采用2020年夏季汛期的地下水位和地下水浓度进行模型拟合。模型长65 m，宽40 m，深度4 m，采用三角形有限单元网格进行离散，在河道区域加密。根据实测地表高程设置模型，垂直分为4层，每层由1069个单元格构成。系统自动生成每天9个模拟时间步长，共计44 d。

图2 FEFLOW中的数字高程模型Fig.2 Digital elevation model in FEFLOW

河岸潜流带中水流采用非稳定流潜水方程描述［26］：

式中：Kxx、Kyy、Kzz是沿坐标轴x、y、z方向的渗透系数（m·d-1），假设坐标轴与渗透系数的方向一致；H为地下水位（m）；h为潜水层厚度（m）；t为时间（d）；ω为单位时间水平面积上的垂向补给量（m3·d-1）；μ为给水度。

溶质反应运移模型采用以下对流弥散方程描述［27］：

式中：c为水中溶质的浓度（mg·L-1）；t为时间（d）;D为水动力弥散系数（m2·s-1)，分为纵向和横向，DL表示纵向弥散系数（m），DT表示横向弥散度（m），可假设其与溶质浓度无关；v为平均地下水流速（m·s-1）；∇为梯度；∇c为浓度梯度；Rc为溶质源汇项（mg·m-3·s-1）。

式中：α为弥散度（m），分为纵向和横向，αL表示纵向弥散度（m），αT表示横向弥散度（m）；D0为分子扩散系数（m2·s-1）。

式中：Rc1为的源汇项；Nlea表示土壤淋滤进入地下水的量（mg·m-3·s-1）；K1表示由各类反应，如反硝化导致浓度降低的衰减系数（d-1）；为浓度；Rc2为的源汇项；Nmin表示有机氮矿化量（mg·m-3·s-1）；K2表示由吸附、硝化等过程导致浓度降低的衰减系数（d-1）；为浓度。

在研究区域内的河流位置设为第一类边界（Dirichlet 边界），根据逐日监测的河水位确定。降水量作为源汇项，降水入渗系数根据区域包气带下垫面条件作为经验值计算。课题组前期研究结果表明［32］，研究河段侧向潜流交换边界位于距离河道10～15 m的位置，该模型前后边界距河流约为20 m，此处河岸地下水与河流之间的交换非常微弱，因此设置为零通量边界。左右边界可能有地下水流入流出，设置为流量边界。垂向潜流交换边界位于河床表面以下1 m 左右［32-33］，该模型下边界位于河床表面以下4 m，基本无垂向交换，因此模型底部设为零通量边界。河水及浓度分别作为溶质模块河流的定浓度边界条件，其余边界的浓度均为0。以2020 年7 月3 日的地下水位和溶质浓度的空间分布作为初始条件，根据监测期内实测水位及溶质浓度对模型进行率定。

运用抽水法测定河床和河岸带沉积物的渗透系数［34］，模型参数［35-39］见表1。

表1 模型参数取值Tab.1 Model parameter values

为了验证模型模拟精度，使用均方根误差（Root Mean Square Error，RMSE）和拟合优度（Coefficient of Determination，R2）作为评价指标，公式如下：

式中：y表示实际测量值；表示实际测量平均值；ym表示模型模拟值；i=1，2，…，n，n表示数据个数。

1.4 模型参数敏感性分析

根据模型模拟污染物溶质特征，选取了5 种对河岸潜流带溶质运移存在影响的参数进行敏感性分析。参数包括：纵向弥散度αL、横向弥散度αT、分子扩散系数D0、一级反应常数K1、K2。局部灵敏度分析是用来分析单个参数变化对数值模型输出结果的影响程度，不考虑参数之间的耦合作用。对一种参数进行灵敏度分析时，保持其他参数不变，将其取值分别增大10%、20%和40%，分别记录不同取值情况下的模拟浓度相对于调试验证后模型结果（基本情景）的变化。

1.5 河岸潜流带去除效果评价

采用Rn反映河岸潜流带内的产生和去除情况［40］：

式中：Chz表示河岸潜流带中浓度表示河水中浓度（mg·L-1）。

基于Kunz等［41］的研究方法，可计算河岸潜流带去除率R。计算公式如下：

2 结果与分析

2.1 水动力学过程

监测期内发生两次较大规模的降雨事件：7 月11—13 日连续3 d 降雨47 mm，8 月8 日降雨64 mm（图3）。受降雨影响，河水位与河岸地下水位在监测期内频繁波动，并且呈现相似的变化趋势，水位变化幅度有所不同。7 月9 日河水位与地下水位小幅升高，监测区内无较大降水，推测可能是由于上游来水导致。在7月11—13日连续降水期间，地下水位达到监测期第一次较大峰值，上涨近0.30 m，河水位略高于地下水位，上涨0.45 m。降水结束后，水位恢复至降雨前7 月10 日的水平。8 月8 日第二次降雨期间，水位较为平稳，仅有小幅度升高。降雨分别导致地下水位及河水位上涨0.20 m、0.27 m。模拟结果表明模型能够较准确得描述地下水位的实际波动过程（图3），模拟和实测的地下水位拟合程度较好，RMSE 为0.04～0.07 m，R2在0.57～0.85 之间。模拟值比实测值相比略有出入，可能是由于不同规模降雨时的河岸入渗能力有所差异，导致地下水位的涨幅稍有不同。

图3 降水量、河水位及河岸监测井实测水位与模拟水位拟合Fig.3 Rainfall，river stage and fitting of measured and simulated water levels in river bank monitoring wells

监测期内频繁的水位波动导致了河水与河岸地下水水力梯度的变化。无降雨时水力梯度为较小的负值，表现为地下水向河流微弱的排泄趋势（图4）。当降雨发生时，河水位上涨幅度较大，导致河水位高于河岸地下水位，水力梯度发生逆转，表现为河水入渗补给地下水，7 月12 日达到监测期内最大值0.047。降雨结束后，河水位下降较快，恢复至降雨前水位，水力梯度为负值，且接近0。

图4 监测期内河水与河岸地下水之间水力梯度随时间变化Fig.4 Change of hydraulic gradient between river water and river bank groundwater during the monitoring period

2.2 氮素浓度变化

图5 不同监测井 浓度的实测值与模拟值拟合Fig.5 Fitting of measured and simulated concentrations in different monitoring wells

图6 不同监测井浓度实测值与模拟值拟合Fig.6 Fitting of measured and simulated concentrations in different monitoring wells

由敏感性分析结果（表2）可知，纵向弥散度αL和横向弥散度αT参数取值增加40%时，引起浓度变化的最大值仅为0.067 mg·L-1，而分子扩散系数D0取值增加40%时，结果变化仍小于0.001 mg·L-1。因此以上3个参数对于及两种溶质的浓度影响较小，参数灵敏度较低。相对而言，的衰减系数K1与的衰减系数K2参数灵敏度较高。当K1增加40%时，峰值浓度降低了2.64 mg·L-1；当K2增加40%时，峰值浓度降低了1.13 mg·L-1。由此可见，反应参数对模型结果影响较大，即河岸带的生物地球化学反应对于氮素归宿具有重要作用。

表2 模型参数敏感性分析结果Tab.2 Parameter sensitivity analysis results of the established model

2.3 相关水环境指标变化

在监测期内，河水pH范围为7.78～8.31，河岸地下水pH 范围为7.52～7.81，均为弱碱性。河水DOC浓度介于8.01～47.3 mg·L-1，河岸监测井DOC浓度介于1.83～29.72 mg·L-1。河岸监测井DOC浓度峰值出现在7 月12 日，DOC 浓度由降水前的11.21～25.18 mg·L-1升至17.27～29.72 mg·L-1（图7）。降雨结束后，DOC浓度呈现下降趋势，在8月2日降至监测期浓度最低值1.37～11.36 mgL-1。河水中DO 充足（7.74±0.34 mg·L-1），可为微生物的好氧呼吸等提供有利条件。河岸地下水中的DO 浓度变化范围为0.2～1.7 mg·L-1，处于缺氧环境，其浓度动态过程与DOC相似，同样在7月12日监测到最大值。

图7 监测期河岸监测井中DOC和DO浓度变化Fig.7 DOC and DO concentration changes in riverbank monitoring wells during the monitoring period

3 讨论

3.1 汛期河岸潜流带中的迁移转化

试验和模拟结果表明，水文过程对河岸潜流带及河流中氮素浓度的变化起了主导作用，但是控制和迁移转化的机制不同。河岸潜流带浓度在连续3 d 强度较小的降雨事件后，小幅上升而后下降，可能是土壤表层中的首先被淋滤出来［46］，在落水期浓度缓慢下降［47］。此外，河岸潜流带DO浓度升高为硝化反应提供了条件（图7），促使了向的转化［48］，导致一定程度河岸潜流带浓度的升高

3.2 汛期河岸潜流带中的迁移转化

8 月8 日的降雨量虽然较大但并未导致水力梯度显著增大，甚至大部分监测井的水力梯度为负，呈现地下水向河水补给，浓度也仅有小幅上升。当降雨结束后，河水与地下水的水力梯度在0附近波动，即河水与地下水交换微弱。缺乏DO 和DOC 进入河岸潜流带的情况不利于矿化作用，浓度也缓慢下降。表明降雨事件通过改变河水-地下水的交换作用，控制营养物质的输入，从而影响河岸潜流带浓度变化。

3.3 牧区河岸带对氮素的源汇效应

夏季汛期降水事件的频发导致了河流及地下水中氮素浓度大幅升高，主要是由于放牧导致大量的氮负荷输入。研究指出在放牧过程中，动物摄入70%～90%的氮会通过粪便和尿液返回牧场［53］。本研究区牧场以牛群为主，而牛群放牧情况下的氮淋失量要高于其他家畜［54］。此外，研究区土壤以砂土为主，因其孔隙较大、渗透性较强，更容易导致淋溶［55］。在整个监测期，河岸潜流带均表现为的产生，而的产生和去除交替出现（图8），因此该处主要考虑的去除效果。根据公式（7）计算得到，7月12日（小雨后）和8月12日（暴雨后）河岸潜流带表现为的产生，是河流氮素的源。小雨后左岸（1L～3L）的产生较右岸（1R～3R）更明显，可能是因为该时段左岸DO 相对右岸较高，不利于的去除（图7）。非降雨期河岸潜流带主要表现为的去除，去除率为14%～86%，是河流氮素的汇，可能因为随着落水期DO 的消耗，河岸潜流带总体趋向厌氧环境（DO＜2 mg·L-1），利于的去除［56］。并且根据方差分析结果（P=0.25＞0.05），河岸潜流带左右岸去除效果较为接近，不具有显著性差异（图8）。

图8 监测期内河岸潜流带产生和去除情况（a）及两岸 去除水平方差分析（b）Fig.8 Production and removal of in riparian hyporheic zones(a)and variance analysis of removal level on both banks(b)during the monitoring period

河岸带植被群落物种较为丰富，可以有效地固定和吸收多余的氮素养分，从而减少了营养物质的流动性，大大增加了营养物质的停留时间，降低淋失的风险［57］。Zuazo 等［58］进行草地植被恢复后与没有植被的裸露土壤相比，发现地下水浓度减少了35%。因此保护河岸带生物多样性，防止水土流失，是牧区氮素消减的有效途径。此外，放牧不仅会改变土壤结构，还会影响微生物过程从而影响氮素的转化和去除，增加淋溶风险［59］，所以加强放牧管理也是控制面源污染的必要措施。

4 结论