连续河水位波动对河床潜流带硝酸盐转化效率的影响

Effect of Continuous River Water Level Fluctuations on Nitrate Conversion Efficiency in Hyporheic Zone

摘要

为了探究连续水位波动下潜流带内硝酸盐的反应迁移规律，构建了包含河床沙丘的垂向二维数值模型.通过考虑不同河水位波动情况、河床坡度以及有氧呼吸、硝化和反硝化过程，系统探讨河床坡度和连续水位波动过程对溶质时空分布以及硝酸盐转化效率的影响.结果表明：河床坡度越大，地表水和周围地下水流之间会发生更快的物质交换，使溶质的浓度变化幅度变小，最终降低潜流带对NO₃ ^‒的转化效率；较高的后续水位峰值会延长地下水潜流路径并使溶质的浓度相对变化幅度越大，然而潜流带对NO₃ ^‒的转化效率会变低；后续水位波动持续的时间会影响溶质的时间响应，但不会影响对NO₃ ^‒的转化效率；不同的后续水位波动延迟时间会影响NO₃ ^‒浓度变化峰值出现的数量，延迟时间越久，越容易出现NO₃ ^‒浓度的多峰现象.

关键词

连续水位波动 / 河床潜流带 / 硝酸盐转化效率 / 反应迁移模拟 / 水文地质.

Key words

continuous water level fluctuations / hyporheic zone in river bed / nitrate conversion efficiency / reactive transport stimulation / hydrogeology

0 引言

当前全球氮污染形式严峻，由于氮肥、农药的大量使用，导致河流、湖泊、土壤、地下水受到大面积的污染，这是人类必须面临的一大全球性的环境问题（Goolsby et al.，2000；Sun et al.，2012；Grant et al.，2018；屈国颖等，2022）.潜流带是地表水和地下水的强烈交互带，是发生生物地球化学反应的活跃场所（Bencala，1983）.潜流带内微生物多种多样的代谢过程以及和非生物物质之间的相互作用影响着潜流带内氮循环，其中反硝化作用是潜流带发挥去除硝酸盐的主要过程（李勇等，2016）.该过程能够促进溶解性无机氮转化为氮气释放，从而离开水流系统，进入大气（Howden et al.，2011）.潜流带作为河流和地下水系统的重要保护屏障，能够有效控制氮污染、净化水质、保护水土环境（吴健等，2006；钱进等，2009；夏继红等，2010）.

连续的河水位波动是现实中的常见现象，大坝等水利设施的调蓄作用会引起下游河流的水位波动更加频繁（Nilsson et al.，2005）.前人研究大多聚焦于动态水位波动过程对潜流交互通量和地下水滞留时间的影响.例如，McCallum and Shanafield（2016）通过构建瞬态水流模型，探究了日变水位波动以及降雨过程对溶质滞留时间分布的影响，发现相比于单峰值水位波动，多峰值水位波动会延长潜流带中地下水的滞留时间.Gomez-Velez et al.（2017）通过构建包含迂曲河岸及单峰值河水位波动的数值模型，发现河岸侧向潜流交互特性主要由含水层水力传导能力和背景地下水流场决定.Harvey et al.（2012）通过分析场地监测数据，发现水位波动峰值越大、持续时间越长，激发潜流过程中溶质和细颗粒物运输的可能性就越大.除此之外，Singh et al.（2020）通过构建二维河床潜流模型探讨了连续水位波动对溶质潜流交换和滞留时间的潜在影响.虽然连续河水位波动对潜流交互的影响已有一定研究，但连续水位波动下潜流过程中的污染物的反应迁移规律尚未被深入探讨.

河水位波动引起的水动力条件的变化，促进了潜流带中的生物地球化学循环，导致潜流带内氮素污染物的衰减更加强烈（Käser et al.，2009；Craig et al.，2010；赵淑凤等，2021）.反硝化反应作为氮循环的核心，是潜流带发挥去除氮素污染物的最重要的途径（Knowles，1982）.前人研究大多聚焦于单峰水位波动过程对潜流带氮循环的影响，Trauth and Fleckenstein（2017）根据实际现场数据建立瞬态地下水流模型，发现河水位波动特征会影响河床的养分供应率，从而在调节潜流带的硝化和反硝化过程中发挥作用.Kessler et al.（2015）建立了河床迁移下潜流带内氮循环的数值模型，发现河床形态迁移对潜流带反硝化作用的影响可以忽略不计.Shuai et al.（2017）通过对多量化参数进行敏感性分析，发现周期性河水位波动特征对潜流带硝酸盐转化率产生了重要的敏感性影响.当前，关于潜流带的硝酸盐转化效率方面已有大量研究，但是尚局限于单峰水位波动过程，对连续水位波动下硝酸盐的迁移转化规律还不清楚，因此仍需进行深入探讨.

本文在前人的研究基础上，通过构建包含河床沙丘的垂向二维数值模型，重点围绕不同的连续河水位波动情况，考虑不同河床坡度以及有氧呼吸、硝化/反硝化过程，系统化地分析了河床坡度和连续水位波动过程对溶质时空分布以及硝酸盐转化效率的影响.本文结果对于更好地评估复杂水位波动和河床形态条件下潜流带对硝酸盐的转化能力至关重要，能为动态河水位下河床含水层污染物去除和河流恢复工作提供重要的参考数据.

1 数值模型的建立

1.1　概念模型

模型区域概化为一个二维的均质各向同性的正弦河床界面（图1a），河床界面压力分布主要受上边界水位变化（H _sup）的影响（如图1b）.正弦河床的振幅Δ为0.1 m，波长λ为1 m，模型长度设置为3个河床周期（L=3λ），深度设置为5个河床周期（D _b=5λ）（Singh et al.， 2019， 2020）.假设河床中所有的溶质都来自于河流的入渗，经过足够长时间后在模型区域内形成稳定的溶质分布.

1.2　水流模型

模型域上边界为正弦河床界面，其在模型中的位置用Z _SWI［L］来进行描述（Singh et al.， 2019）：

Z_{\mathrm{S}\mathrm{W}\mathrm{I}} = \frac{∆}{\mathrm{2}} \mathrm{s} \mathrm{i} \mathrm{n} \left(\frac{\mathrm{2}\mathrm{\pi }x}{\lambda }\right),（1）

其中x［L］为水平坐标，y［L］为垂直坐标.Δ［L］为河床振幅，λ［L］为河床波长.假设河床含水层为均质各向同性，且水力参数不会发生变化.

模型中地下水运动方程使用水均衡方程来描述：

\frac{\partial \left(\rho \theta \right)}{\partial t} + \nabla \bullet \left(\rho \mathit{q}\right) = \mathrm{0},（2）

\mathit{q} = - \frac{K}{\rho \mathrm{g}} \bullet ( \nabla p + \rho \mathrm{g} \nabla z ),（3）

其中z［L］代表位置高程，p［MT^-2L^-1］代表的是压力水头， q ［LT^-1］代表达西流速场矢量，K［ML^-1］为渗透系数，θ［-］为孔隙度，ρ［ML^-3］为流体密度.

河床表面的水头分布是关于大尺度和小尺度河床地形的线性组合（Wörman et al.， 2006； Stonedahl et al.，2010）.使用以下方程来刻画河床表面水头的时空变化：

h_{\mathrm{S}\mathrm{W}\mathrm{I}} \left(x,t\right) = - S x + \left[H_{\mathrm{s}\mathrm{u}\mathrm{p}}\left(t\right)-Z_{\mathrm{S}\mathrm{W}\mathrm{I}}\left(x\right)\right] + \frac{\mathrm{2}{h}_{\mathrm{d}}\left(t\right)}{\mathrm{\Delta }} Z_{\mathrm{S}\mathrm{W}\mathrm{I}} ( x + \frac{\lambda }{\mathrm{4}} ),（4）

其中S［-］为河道坡度，H _sup（t）［L］为河水位波动过程，h _d（t）［L］为动态水位波动强度（Elliott and Brooks，1997）：

h_{\mathrm{d}} \left(t\right) = \mathrm{0.28} \frac{U_{\mathrm{S}}{(t)}^{\mathrm{2}}}{\mathrm{2}\mathrm{g}} \left\{\begin{array}{c}{\left(\frac{∆}{\mathrm{0.34}{H}_{\mathrm{s}\mathrm{u}\mathrm{p}}(t)}\right)}^{\mathrm{3}/\mathrm{8}},\frac{∆}{H_{\mathrm{s}\mathrm{u}\mathrm{p}}(t)}\le \mathrm{0.34}\\ {\left(\frac{∆}{\mathrm{0.34}{H}_{\mathrm{s}\mathrm{u}\mathrm{p}}(t)}\right)}^{\mathrm{3}/\mathrm{2}},\frac{∆}{H_{\mathrm{s}\mathrm{u}\mathrm{p}}(t)}>\mathrm{0.34}\end{array}\right.,（5）

其中，河流平均流速由Chezy方程定义 : U_{\mathrm{s}} \left(t\right) = M^{-\mathrm{1}} H_{\mathrm{s}} {(t)}^{\mathrm{2}/\mathrm{3}} S^{\mathrm{1}/\mathrm{2}}，M［L^-1/3T］为曼宁系数.使用H _generic（t）［L］来描述两个水位过程重叠部分的水位波动：

,（6）

其中H ₀［L］为初始基流条件下的河水位，H _p［L］为河流波动的最大水位变化，t _p［T］为河流从初始水位到达峰值水位的时间，t _lag［T］为后续水位波动过程的发生时间，t _d［T］为水位波动过程的持续时间， \omega = \mathrm{2} \mathrm{\pi } / t_{\mathrm{p}} [ \mathrm{T} - \mathrm{1} ]为水位波动发生的频率， \delta = \omega \mathrm{c} \mathrm{o} \mathrm{t} ( \omega t_{\mathrm{p}} / \mathrm{2} ) [ \mathrm{T} - \mathrm{1} ]为水位峰值偏度.使用H _sup（t）来刻画连续水位波动过程：

H_{\mathrm{s}\mathrm{u}\mathrm{p}} \left(t\right) = H_{\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}} \left(t,H_{\mathrm{p}\mathrm{1}},t_{\mathrm{p}\mathrm{1}},t_{\mathrm{d}\mathrm{1}},\mathrm{0}\right) + H_{\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}} \left(t,H_{\mathrm{p}\mathrm{2}},t_{\mathrm{p}\mathrm{2}},t_{\mathrm{d}\mathrm{2}},t_{\mathrm{l}\mathrm{a}\mathrm{g}}\right) - H_{\mathrm{0}},（7）

其中，H _p1、t _p1和t _d1描述了前一个发生的水位波动过程，H _p2、t _p2、t _d2和t _lag描述了后续发生的水位波动过程.河床形态沿河流周期性重复，设置左右边界的周期性边界条件为：

p（- L，y，t）=p（2L，y，t）+Δp，（8）

底部边界为隔水边界：

\mathit{n} \bullet \rho \mathit{q} = \mathrm{0},（9）

其中 n 为外法线方向.

1.3　多组分反应迁移模型

水位波动将多组分溶质带入潜流带，促进潜流带中生物地球化学反应的发生（Williams et al.，2010）.本文考虑了三种主要的化学反应：有氧呼吸、硝化和反硝化过程.具体的化学反应式如下：

\mathrm{C} {\mathrm{H}}_{\mathrm{2}} \mathrm{O} + {\mathrm{O}}_{\mathrm{2}} \to \mathrm{C} {\mathrm{O}}_{\mathrm{2}} + {\mathrm{H}}_{\mathrm{2}} \mathrm{O},（10）

\mathrm{5} \mathrm{C} {\mathrm{H}}_{\mathrm{2}} \mathrm{O} + \mathrm{4} \mathrm{N} {\mathrm{O}}_{\mathrm{3}}^{-} + \mathrm{4} {\mathrm{H}}^{+} \to \mathrm{5} \mathrm{C} {\mathrm{O}}_{\mathrm{2}} + \mathrm{2} {\mathrm{N}}_{\mathrm{2}} + \mathrm{7} {\mathrm{H}}_{\mathrm{2}} \mathrm{O},（11）

\mathrm{N} {\mathrm{H}}_{\mathrm{4}}^{+} + \mathrm{2} {\mathrm{O}}_{\mathrm{2}} \to \mathrm{N} {\mathrm{O}}_{\mathrm{3}}^{-} + \mathrm{2} {\mathrm{H}}^{+} + {\mathrm{H}}_{\mathrm{2}} \mathrm{O},（12）

其中河水作为DOC（溶解性有机碳，Dissolved Organic Carbon）、DO（溶解氧，Dissolved Oxygen）、硝酸根离子（NO₃ ^-）和铵根离子（NH₄ ⁺）四种主要反应溶质的输入源，通过河床表面和潜流过程进入模型区域.其中DOC的微生物降解被视为主要反应（Hunter et al.，1998），遵循一级反应动力学（Bardini et al.，2012）：

r_{\mathrm{D}\mathrm{O}\mathrm{C}} = k_{\mathrm{D}\mathrm{O}\mathrm{C}} \bullet \left[\mathrm{D}\mathrm{O}\mathrm{C}\right],（13）

其中k _DOC为一级反应速率常数，［DOC］表示DOC的实际活度，活度系数设为1.硝化反应用常规的二级双分子反应动力学描述：

r_{\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{+}} = k_{\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{+}} \bullet [ \mathrm{N} {\mathrm{H}}_{\mathrm{4}}^{+} ] \bullet \left[{\mathrm{O}}_{\mathrm{2}}\right],（14）

其中k _NH4 ⁺ 为二阶反应速率常数，［NH₄ ⁺］为NH₄ ⁺的实际活度，活度系数设为1.假设有氧呼吸和硝化作用同时发生，而反硝化作用只在氧浓度低于某一限值时才会发生.当DOC降解时，释放的电子转移到O₂和NO₃ ^-终端电子受体，可以计算出第i个电子受体的还原率（i=1是O₂，i=2是NO₃ ^-）（Zheng and Cardenas，2018）：

r_{\mathrm{r}\mathrm{e}\mathrm{d},i} = f_i \bullet r_{\mathrm{D}\mathrm{O}\mathrm{C}} \bullet {\beta }_i,（15）

其中β _i是每摩尔氧化DOC转移电子的摩尔数与每摩尔还原DOC转移电子的摩尔数之比，O₂转移的电子比为1，NO₃ ^-转移的电子比为0.8.f_i 是指还原反应消耗电子的比例，用于终端电子受体，用简化的Monod公式定义（Hunter et al.， 1998； Bardini et al.， 2012）：

f_i = \left(\mathrm{1}-{\sum }_{n=\mathrm{0}}^{i-\mathrm{1}}{f}_n\right) \bullet {\alpha }_i,（16）

,（17）

式中C_i 和C_{i ，lim}分别为第i个电子受体的摩尔浓度和摩尔极限浓度，f ₀=0.α_i 是一个无因次参数，它假设每个电子受体有一个极限浓度（C_{i ，lim}）.当该电子受体（C_i ）的浓度超过阈值时，相应的反应速率与C_i 无关.否则，速率与C_i 成正比.有些物质在反应中起双重作用，既作为反应物又作为生成物，因此定义这四种反应物的净反应速率为：

R_{\mathrm{D}\mathrm{O}\mathrm{C}} = - r_{\mathrm{D}\mathrm{O}\mathrm{C}},（18）

R_{{\mathrm{O}}_{\mathrm{2}}} = - r_{\mathrm{r}\mathrm{e}\mathrm{d},\mathrm{1}} - \mathrm{2} r_{\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{+}},（19）

R_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}^{-}} = - r_{\mathrm{r}\mathrm{e}\mathrm{d},\mathrm{2}} + r_{\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{+}},（20）

R_{\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{+}} = - r_{\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{+}},（21）

负号表示反应物被消耗，正号表示反应物的生成.含水层中溶质的反应迁移过程使用对流弥散方程（Advection-Dispersion Equation，ADE）来计算：

\theta \frac{\partial C_{\mathrm{i}}}{\partial t} = \nabla \bullet ( \mathit{D} \nabla C_{\mathrm{i}} ) - \nabla ( \mathit{q} C_{\mathrm{i}} ) + R_i ( i = \mathrm{D} \mathrm{O} \mathrm{C} , \mathrm{D} \mathrm{O} , \mathrm{N} {\mathrm{H}}_{\mathrm{4}}^{+} , \mathrm{N} {\mathrm{O}}_{\mathrm{3}}^{-} ),（22）

其中C_i ［ML^-3］为不同溶质的摩尔浓度，R_i ［ML^-3T^-1］为不同溶质的反应速率， D =｛D_ij ｝为水动力弥散系数张量：

D_{ij} = {\alpha }_{\mathrm{T}} \left|\mathit{q}\right| {\delta }_{\mathrm{i}\mathrm{j}} + \left({\alpha }_{\mathrm{L}}-{\alpha }_{\mathrm{T}}\right) \frac{q_i{q}_j}{\left|\mathit{q}\right|} + \frac{\theta }{{\xi }_{\mathrm{m}}} D_{\mathrm{m}} （23）

其中α _T［L］和α _L［L］为横向和纵向弥散度，D _m［L²T^-1］为有效分子扩散系数，ξ _m=θ ^-1/3［-］为孔隙曲折度，δ_ij ［-］为克罗拉格函数（Kronecker function）.河流为各反应组分的来源，因此河床边界为开放边界：

,（24）

\frac{\partial C_{\mathrm{i}}\left(x,y,t\right)}{\partial \mathit{n}} = \mathrm{0} \left(i=\mathrm{D}\mathrm{O}\mathrm{C},\mathrm{D}\mathrm{O},\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{+},\mathrm{N}{\mathrm{O}}_{\mathrm{3}}^{-}\right) \mathit{n} \bullet \mathit{q} \ge \mathrm{0},（25）

其中 n 为外法线方向.当河床为流入边界时，其为溶质的定浓度边界（式24）；当河床为流出边界时，其为零通量边界（式25）.为了避免边界效应对结果的影响，且模拟无限长的河床含水层，左右边界处溶质为周期性边界条件为：

C_i \left(-L,y,t\right) = C_i \left(\mathrm{2}L,y,t\right),（26）

底部的无流动边界条件为：

\mathit{n} \bullet ( \mathit{D} \nabla C_i - \mathit{q} C_i ) = \mathrm{0},（27）

其中 n 为外法线方向.

1.4　不同连续水位波动情景和参数设置

本文构建两个连续水位波动对潜流流场和溶质循环的影响，考虑不同波动延迟、峰值和持续时间三种影响因素下的连续水位波动情景：①两个水位波动过程发生的延迟时间不同（t _lag = t_{\mathrm{p}\mathrm{1}} + \eta ( t_{\mathrm{d}\mathrm{1}} - t_{\mathrm{p}\mathrm{1}} )，其中η为前后两个水位波动过程启动时间间隔）；②水位波动的峰值不同（H _p1-H _p2）；③水位波动的持续时间不同（t _d1 -t _d2），如图2所示.图2中情景a为基准情景，当η=1/4、2/4和3/4时两个水位过程发生叠加，当η=4/4时，两个水位过程相互不叠加也不干涉.情景b和c中同样考虑不同后续波动延迟时间，同时情景b中后续水位波动峰值相对于情景a更高，但持续时间相同.情景c中后续水位波动持续时间比情景a更长，但波动峰值相同.情景b和情景c可以与情景a进行对比分析.

此外，基于每个水位波动情景还考虑了不同河床坡度S分别为0.001、0.01、0.05的情况，用来综合分析水位波动和河床形态的影响.综上所述，本文在三组模拟情景下共构建36个数值模型，综合对比分析了不同后续水位波动和河床坡度对结果的影响.

本文通过对不同波动过程分别进行溶质示踪，来衡量单独的水位波动过程引起的潜流交互对溶质迁移转化的贡献，还可以探讨水位叠加时段的影响.如图3所示，C ₁表示只对前一个水位波动的时段进行溶质示踪，C ₂表示对两个水位波动的叠加时段进行示踪，C ₃表示只对后续水位波动时段进行示踪.

本文用流量乘以流入和流出潜流带的NO₃ ^-浓度，来刻画进入和离开体系的NO₃ ^-质量.定义进入潜流带的NO₃ ^-质量为：

M_{\mathrm{i}\mathrm{n}} = {\int }_a^b \left[{\int }_{\partial {\mathrm{\Omega }}_{\mathrm{I}\mathrm{N}}}\left(\mathit{n}\bullet \mathit{q}\right)C_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}^{-}}\mathrm{d}x\right] \mathrm{d} t，（28）

式中∂Ω_IN是河床表面的流入边界.定义离开潜流带的NO₃ ^-总质量为：

m_{\mathrm{o}\mathrm{u}\mathrm{t}} = {\int }_{\partial {\mathrm{\Omega }}_{\mathrm{O}\mathrm{U}\mathrm{T}}} \left(\mathit{n}\bullet \mathit{q}\right) C_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}^{-}} \mathrm{d} x，（29）

式中∂Ω_OUT是河床表面的流出边界.通过绘制NO₃ ^-的穿透曲线，可以分析不同水位波动下潜流带对NO₃ ^-转化效率的影响.定义NO₃ ^-的穿透曲线（BTC）为：

\mathrm{B} \mathrm{T} \mathrm{C} \left(t\right) = \frac{{\int }_{\mathrm{0}}^t{m}_{\mathrm{o}\mathrm{u}\mathrm{t}}\mathrm{d}t}{M_{\mathrm{i}\mathrm{n}}} .（30）

本文通过多物理场仿真模拟软件COMSOL Multiphysics来构建数值模型.使用自由三角形网格来划分整个模型区域.为了在保证正确性的基础上减少计算量，考虑到溶质的迁移过程只发生在模型域上半部分，因此着重细化了河床上半部的网格，并在河床边界处划分边界层网格进一步提高模型精度.整个模型域共划分成213 769个三角形网格，111 245个网络顶点，1 039个边单元.模型的计算步骤分为两步：①先使用稳态求解器求解出稳定状态下的流场和溶质场；②使用瞬态求解器将步骤1得到的稳态结果作为瞬态模拟的初始值，得到瞬态的流场和溶质场.表1为模型中的各参数符号及默认取值.

2 结果分析

2.1　流场及潜流带范围

本文在水位波动过程中共选取了如图4所示6个时刻点对流场和溶质场进行分析，其中：a为t=0时初始基流H ₀条件下水位，b为前一个水位波动达到峰值时的水位；从c（t=t _lag）开始，两个水位波动开始叠加；d为后一水位波动达到峰值时的水位；e和f分别表示前后的水位波动过程结束时的水位.

图5为情景a下η=3/4时不同时间点流场的时空分布图，颜色图例表示的是达西流速的大小.黑线代表潜流带的范围.图5a为t=0时初始基流条件下的流场，从t=t _lag（图5c）开始，流场受到两个水位波动叠加的影响，河床内流线的渗透深度会逐渐扩大，直至最后恢复到初始状态.由于前一个水位波动过程的峰值更高，当t=t _p1时流线的渗透深度（图5b）要比t=t _p2时后续水位波动达到峰值（图5d）的渗透深度更大.由此可以看出，水位波动的峰值大小和持续时间决定流场的动态变化.当两个连续水位波动过程发生叠加，会促进水流的入渗，且入渗深度受水位峰值的影响.

本文使用Triska et al.（1989）对潜流带范围的量化方法，即将含水层中50%的孔隙水来自地表水的区域定义为潜流带.结果表明，潜流带的范围随时间呈现不同程度的波动起伏.这种波动变化主要出现在t=t _p1（图5b）和t=t _p2（图5d）之间，都与连续水位波动的两次水位上升有关.除此之外，流场中存在达西流速为零或接近零的区域，该区域一般被称作停滞点区域（图5中黄色点域）.图5中停滞点区域很大程度上由河水位波动决定，随着流线渗透深度的增大停滞点逐渐向含水层深处运动.停滞区作为潜流带发生生物地球化学反应的“热点”区域，能够创造合适的氧化还原条件，促进硝化和反硝化反应的发生.受垂直于河床的水力梯度的影响，流场和动态停滞区的动态变化会很大程度上改变潜流带内的硝化和反硝化过程，对NO₃ ^-的迁移转化产生重要影响.

2.2　净反硝化速率及NO₃ ^-的时空分布

在连续水位波动过程中，河床内反应速率及浓度分布的时空变化较大.图6为坡度S=0.01以及延迟时间η=3/4时，三种情景下不同时间点NO₃ ^-（左三列）及净反硝化速率（右三列）的时空分布.随着水位的上升和下降，不同情景下NO₃ ^-以及净反硝化速率的时空分布呈现出相似的扩大和收缩趋势.当t=t _p1（图6b）时，前一个水位波动过程峰值达到最高，NO₃ ^-的穿透深度以及净反硝化速率达到最大.从t=t _lag（图6c）开始，两个水位波动过程开始发生叠加，溶质会继续向深处渗透，但渗透深度受后续水位波动的影响.情景b的后续水位波动峰值更高，因此较其他两个情景溶质渗透深度更深，反应速率更强.情景c后续水位波动的持续时间更长，溶质到达峰值所需的时间就更长.情景a和情景b下溶质到达峰值所需时间为1.06 d，而情景c所需的时间为1.13 d.

水位波动过程会驱动不同深度的径流路径，导致溶质浓度的梯度变化.随着水位上升溶质渗入地下水，高浓度反应物在靠近河床表面的区域发生耗氧反应.随着溶质渗透深度增大，来自河流的NO₃ ^-以及硝化反应产生的NO₃ ^-在含水层深处发生反硝化反应.反硝化区域在硝化区域周围形成带状分布，并在河床中部逐渐变浅，呈月牙状分布.这些结论表明，在不同水位波动影响因素下，不同深度的径流路径上硝酸盐的迁移转化特征会有很大不同，需要进行更细致的探讨.因此，后续将对含水层的硝酸盐转化效率进行定量分析.

2.3　不同情景下反应物浓度的相对变化

图7为情景a下，不同的连续水位波动过程以及不同的河床坡度下反应物浓度与基流条件（H=H ₀）的相对变化；图8和图9分别为情景b及情景c的结果图.横向每一行代表不同坡度的结果，纵向每一列代表后续水位波动不同延迟时间的结果，灰线代表的是水位波动过程.在所有情景下，河床内反应物浓度的相对变化都大于0，说明与基流条件相比在水位波动条件下河床内溶质浓度更大.各反应物浓度与水位波动过程的变化趋势基本相一致，其中NH₄ ⁺浓度变化幅度最小.DOC和NO₃ ^-的浓度变化峰值受坡度的影响，随坡度的增大而减小.当坡度很小时（S=0.001），溶质通过浅表较短的径流路径，反应物会快速得到补给，故浓度变化幅度较大.随着坡度增大，溶质逐渐渗透到河床深处，氧化反应和硝化反应消耗大量的DOC和NH₄ ⁺，此时浓度变化峰值较低.在河床深部发生的反硝化反应促进NO₃ ^-转化为N₂释放，NO₃ ^-浓度逐渐降低.

不同水位波动特征会对河床中的溶质浓度产生较大的影响.情景b（图8）与情景a（图7）相比后续水位峰值更高，各溶质的浓度相对变化幅度更大.情景c（图9）与情景a（图7）相比后续水位波动持续时间更长，各溶质浓度的相对变化波动更长.此外，在所有情景中，后续水位波动不同的延迟时间（η）会影响NO₃ ^-浓度变化峰值出现的数量，随着延迟时间的增大，越容易出现多峰.当坡度为0.01时最明显，在η=1/4时，NO₃ ^-浓度变化为双峰，随着延迟时间的增大，浓度变化逐渐出现三峰甚至多峰.

图7~9中NO₃ ^-浓度变化出现多峰的现象，受河床内发生的硝化和反硝化反应的影响.在初始时段，河床表层的DOC和NH₄ ⁺会与DO反应生成大量的NO₃ ^-.在河床深层，NO₃ ^-发生反硝化反应被消耗.由于初始时段氧气充足，硝化反应占主导，因此表现出大量NO₃ ^-的生成，这是NO₃ ^-浓度出现第一个峰的原因.此时DO被大量消耗，反硝化区逐渐上升且范围逐渐扩大，反硝化反应占主导作用.当前一个水位波动逐渐消退，后续水位过程还未开始，潜流层恢复到基流条件，此时硝化和反硝化反应处于互相竞争的状态.随着流速的减小，少量溶质进入潜流带，NO₃ ^-浓度会出现第二次和第三次的上升，但此时氧气会迅速被消耗，反硝化再次占主导，消耗生成的NO₃ ^-.当后续水位波动开始，溶质大量渗入河床，硝化反应再次占主导，NO₃ ^-浓度变化会出现第四次上升.

结果表明，河床内反应物浓度的时空分布高度依赖于河床地形和周围地下水流量引起的局部压力的变化.坡度影响浓度变化峰值的量级，随着坡度值的增加，地表水和周围地下水流之间发生更快的物质交换，促进溶质不断渗入地下水.水位波动特征可以很大程度上决定河床中溶质的反应迁移，后续水位峰值越高，各溶质的浓度相对变化幅度越大.后续水位波动持续时间越长，各溶质的浓度相对变化波长越长.值得关注的是，延迟时间会影响NO₃ ^-浓度变化峰值出现的数量，延迟时间越大，越容易出现多峰.

2.4　不同情景和坡度对NO₃ ^‒转化效率的影响

通过溶质示踪，可以深入了解NO₃ ^-在河床中的迁移和转化.图10、图11和图12分别为情景a、情景b及情景c下NO₃ ^-的穿透曲线示意图.横向每一行代表不同坡度的结果，纵向每一列代表后续水位波动不同延迟时间的结果.

在所有情景下，C ₁、C ₂、C ₃示踪的NO₃ ^-穿透曲线值大致呈现逐渐递减的趋势，说明连续的水位波动过程会逐渐促进潜流带对NO₃ ^-的转化效率.坡度越大，穿透曲线越接近于1，说明NO₃ ^-的转化效率越小.原因在于坡度越大，河床横向边界的流速越快，增加了潜流带对NO₃ ^-的转化负担，故对NO₃ ^-的转化效率会变小.除此之外，随着延迟时间η值的增大，不同示踪段的NO₃ ^-穿透曲线值相差不超过0.1，说明对NO₃ ^-转化效率影响在10%左右.在前一个水位波动下（C ₁线）穿透曲线值随着η值的增大而减小，说明此时潜流带的转化效率会逐渐增大.以情景a为例（图10），坡度为0.01时（第二行），NO₃ ^-的转化效率由32%增长为35%、38%和41%；叠加水位影响下（C ₂线）对NO₃ ^-的转化效率影响很小，基本没有变化；后续水位波动下（C ₃线）穿透曲线值随着η值的增大而增大，说明此时潜流带的转化效率会逐渐减小.同样以情景a为例（图10），坡度为0.01时（第二行），NO₃ ^-的转化效率由68%减小为67%和62%.原因在于C ₁、C ₂、C ₃示踪的是一个连续的水文过程，延迟时间会影响不同示踪段溶质的传质速率，从而对NO₃ ^-的转化效率产生影响.

不同水位波动情景会影响潜流带内NO₃ ^-的转化效率.情景b（图11）与情景a（图10）相比后续水位峰值更高，两个水位波动过程叠加影响下（C ₂线），随着穿透曲线的值逐渐增大，潜流带对NO₃ ^-的转化效率逐渐降低.原因在于连续的高峰值水位波动会增强上边界的水头压力，生成较长的地下水径流路径，导致潜流带对NO₃ ^-的转化效率变低.情景c（图12）与情景a（图10）相比后续水位波动持续时间更长，后续水位波动影响下（C ₃线），穿透曲线逐渐趋于平缓，但穿透曲线值相差不大，说明对NO₃ ^-的转化效率影响很小，两种情况下NO₃ ^-穿透曲线值相差不大.原因在于水位波动持续时间越长，溶质的时间响应会越长，穿透曲线越会趋于平缓，但不会影响潜流带对NO₃ ^-的转化效率.

结果表明，不同的水位波动特征和河床地形会影响潜流带对NO₃ ^-的转化效率.坡度越大，河床横向边界的流速越快，会降低NO₃ ^-的转化效率.延迟时间会影响不同示踪段溶质的传质速率，从而对NO₃ ^-的转化效率产生影响，影响范围基本保持在10%以内.后续水位峰值越高，潜流带对NO₃ ^-的转化效率会降低.后续水位波动持续时间只影响溶质的时间响应，不会影响潜流带对NO₃ ^-的转化效率.

3 结论

本文围绕潜流带内发生的有氧呼吸、硝化和反硝化过程，探讨了河床坡度和不同连续水位波动对潜流交互过程中溶质的时空分布以及NO₃ ^-转化效率的影响.结果表明：

（1）坡度越大，地表水和周围地下水流之间发生更快的物质交换，促进溶质的入渗，从而影响溶质浓度变化峰值的量级，增加了潜流带对NO₃ ^-的转化负担.

（2）后续水位峰值越高，会增强上边界的水头压力，形成较长的地下水径流路径，溶质的浓度相对变化幅度越大，潜流带对NO₃ ^-的转化效率变低.

（3）后续水位波动持续时间越长，各溶质的浓度相对变化波长越长，只影响溶质的时间响应，不会影响潜流带对NO₃ ^-的转化效率.

（4）后续水位波动延迟时间越大，NO₃ ^-浓度变化越容易出现多峰.延迟时间还会影响不同示踪段溶质的传质速率，从而影响NO₃ ^-的转化效率.

本文的研究结果提高了在复杂水位波动和河床形态条件下对潜流带NO₃ ^-转化能力的认识，有助于预测连续水文情景对河流水生态健康的影响，为洪水过后河流恢复工作提供重要的参考数据，可有效地开展对河流生态恢复的管理工作，制定适宜的生态修复策略.

参考文献

[1] Bardini, L., Boano, F., Cardenas, M. B., et al., 2012. Nutrient Cycling in Bedform Induced Hyporheic Zones. Geochimica et Cosmochimica Acta, 84: 47-61. https://doi.org/10.1016/j.gca.2012.01.025

[2] Bencala, K. E., 1983. Simulation of Solute Transport in a Mountain Pool-and-Riffle Stream with a Kinetic Mass Transfer Model for Sorption. Water Resources Research, 19(3): 732-738. https://doi.org/10.1029/WR019i003p00732

[3] Craig, L.S., Bahr, J. M., Roden, E. E., 2010. Localized Zones of Denitrification in a Floodplain Aquifer in Southern Wisconsin, USA. Hydrogeology Journal, 18(8): 1867-1879. https://doi.org/10.1007/s10040-010-0665-2

[4] Elliott, A. H., Brooks, N. H., 1997. Transfer of Nonsorbing Solutes to a Streambed with Bed Forms: Theory. Water Resources Research, 33(1): 123-136. https://doi.org/10.1029/96wr02784

[5] Gomez-Velez, J. D., Wilson, J. L., Cardenas, M. B., et al., 2017. Flow and Residence Times of Dynamic River Bank Storage and Sinuosity-Driven Hyporheic Exchange. Water Resources Research, 53(10): 8572-8595. https://doi.org/10.1002/2017wr021362

[6] Goolsby, D. A., Battaglin, W. A., Aulenbach, B. T., et al., 2000. Nitrogen Flux and Sources in the Mississippi River Basin. Science of the Total Environment, 248(2-3): 75-86. https://doi.org/10.1016/s0048-9697(99)00532-x

[7] Grant, S. B., Azizian, M., Cook, P., et al., 2018. Factoring Stream Turbulence into Global Assessments of Nitrogen Pollution. Science, 359(6381): 1266-1269. https://doi.org/10.1126/science.aap8074

[8] Harvey, J. W., Drummond, J. D., Martin, R. L., et al., 2012. Hydrogeomorphology of the Hyporheic Zone: Stream Solute and Fine Particle Interactions with a Dynamic Streambed. Journal of Geophysical Research: Biogeosciences, 117(G4): G00N11. https://doi.org/10.1029/2012jg002043

[9] Howden, N. J. K., Burt, T. P., Worrall, F., et al., 2011. Nitrate Pollution in Intensively Farmed Regions: What are the Prospects for Sustaining High-Quality Groundwater? Water Resources Research, 47(6): W00L02. https://doi.org/10.1029/2011wr010843

[10] Hunter, K. S., Wang, Y. F., Van Cappellen, P., 1998. Kinetic Modeling of Microbially-Driven Redox Chemistry of Subsurface Environments: Coupling Transport, Microbial Metabolism and Geochemistry. Journal of Hydrology, 209(1-4): 53-80. https://doi.org/10.1016/s0022-1694(98)00157-7

[11] Käser, D. H., Binley, A., Heathwaite, A. L., et al., 2009. Spatio-Temporal Variations of Hyporheic Flow in a Riffle-Step-Pool Sequence. Hydrological Processes, 23(15): 2138-2149. https://doi.org/10.1002/hyp.7317

[12] Kessler, A. J., Cardenas, M. B., Cook, P. L. M., 2015. The Negligible Effect of Bed Form Migration on Denitrification in Hyporheic Zones of Permeable Sediments. Journal of Geophysical Research: Biogeosciences, 120(3): 538-548. https://doi.org/10.1002/2014jg002852

[13] Knowles, R., 1982. Denitrification. Microbiological Reviews, 46(1): 43-70. https://doi.org/10.1128/mr.46.1.43-70.1982

[14] Li, Y., Zhang, W.W., Yuan, J.H., et al., 2016. Research Advances in Flow Patterns and Nitrogen Transformation in Hyporheic Zones. Journal of Hohai University (Natural Sciences), 44(1): 1-7 (in Chinese with English abstract).

[15] McCallum, J. L., Shanafield, M., 2016. Residence Times of Stream-Groundwater Exchanges Due to Transient Stream Stage Fluctuations. Water Resources Research, 52(3): 2059-2073. https://doi.org/10.1002/2015wr017441

[16] Nilsson, C., Reidy, C. A., Dynesius, M., et al., 2005. Fragmentation and Flow Regulation of the World’s Large River Systems. Science, 308(5720): 405-408. https://doi.org/10.1126/science.1107887

[17] Qian, J., Wang, C., Wang, P.F., et al., 2009. Research Progresses in Purification Mechanism and Fitting Width of Riparian Buffer Strip. Advances in Water Science, 20(1): 139-144 (in Chinese with English abstract).

[18] Qu, G.Y., Li, M.J., Zheng, J.H., et al., 2022. The Promoting Effect and Mechanism of Nitrogen Conversion in the Sediments of Polluted Lake on the Degradation of Organic Pollutants. Earth Science, 47(2): 652-661 (in Chinese with English abstract).

[19] Shuai, P., Cardenas, M. B., Knappett, P. S. K., et al., 2017. Denitrification in the Banks of Fluctuating Rivers: The Effects of River Stage Amplitude, Sediment Hydraulic Conductivity and Dispersivity, and Ambient Groundwater Flow. Water Resources Research, 53(9): 7951-7967. https://doi.org/10.1002/2017wr020610

[20] Singh, T., Gomez-Velez, J. D., Wu, L. W., et al., 2020. Effects of Successive Peak Flow Events on Hyporheic Exchange and Residence Times. Water Resources Research, 56(8): e2020WR027113. https://doi.org/10.1029/2020wr027113

[21] Singh, T., Wu, L. W., Gomez-Velez, J. D., et al., 2019. Dynamic Hyporheic Zones: Exploring the Role of Peak Flow Events on Bedform-Induced Hyporheic Exchange. Water Resources Research, 55(1): 218-235. https://doi.org/10.1029/2018WR022993

[22] Stonedahl, S. H., Harvey, J. W., Wörman, A., et al., 2010. A Multiscale Model for Integrating Hyporheic Exchange from Ripples to Meanders. Water Resources Research, 46(12): W12539. https://doi.org/10.1029/2009wr008865

[23] Sun, B., Zhang, L.X., Yang, L.Z., et al., 2012. Agricultural Non-Point Source Pollution in China: Causes and Mitigation Measures. Ambio, 41(4): 370-379. https://doi.org/10.1007/s13280-012-0249-6

[24] Trauth, N., Fleckenstein, J. H., 2017. Single Discharge Events Increase Reactive Efficiency of the Hyporheic Zone. Water Resources Research, 53(1): 779-798. https://doi.org/10.1002/2016WR019488

[25] Triska, F. J., Kennedy, V. C., Avanzino, R. J., et al., 1989. Retention and Transport of Nutrients in a Third-Order Stream in Northwestern California: Hyporheic Processes. Ecology, 70(6): 1893-1905. https://doi.org/10.2307/1938120

[26] Williams, D. D., Febria, C. M., Wong, J. C. Y., 2010. Ecotonal and Other Properties of the Hyporheic Zone. Fundamental and Applied Limnology, 176(4): 349-364. https://doi.org/10.1127/1863-9135/2010/0176-0349

[27] Wörman, A., Packman, A. I., Marklund, L., et al., 2006. Exact Three-Dimensional Spectral Solution to Surface-Groundwater Interactions with Arbitrary Surface Topography. Geophysical Research Letters, 33: L0740. https://doi.org/10.1029/2006gl025747

[28] Wu, J., Huang, S.F., Tang, H., et al., 2006. Review of Research on Ecosystem Health in Riverine Phreatic Zones. Water Resources Protection, 22(5): 5-8, 27 (in Chinese with English abstract).

[29] Xia, J.H., Lin, J.Q., Yao, L., et al., 2010. Edge Structure and Edge Effect of Riparian Zones. Journal of Hohai University (Natural Sciences), 38(2): 215-219 (in Chinese with English abstract).

[30] Zhao, S.F., Liu, H., Zhao, L., et al., 2021. Responses of Different Iron and Nitrogen Transformation Functional Microorganisms to Fe(Ⅱ) Chemical Oxidation. Earth Science, 46(4): 1481-1489 (in Chinese with English abstract).

[31] Zheng, L. Z., Cardenas, M.B., 2018. Diel Stream Temperature Effects on Nitrogen Cycling in Hyporheic Zones. Journal of Geophysical Research: Biogeosciences, 123(9): 2743-2760. https://doi.org/10.1029/2018jg004412

基金资助

国家自然科学基金面上项目(42272290)

国家自然科学基金重点项目(41830862)

国家自然科学基金项目(U23A2042)