砂质盆地多级次地下水流系统中硝酸盐水文遗留时间模拟及管控启示

Simulation of Nitrate Hydrologic Residence Time in a Hierarchical Groundwater Flow System of Sandy Basin and Implications from Management

摘要

以小型砂质盆地为例，基于二维Tóth模型探究了硝酸盐在多级次地下水流系统中的水文遗留时间特征.随着入渗强度增大，水流系统从单一区域系统、嵌套系统向单一局部水流系统转化，盆地硝酸盐平均水文遗留时间递减，但较长遗留时间（250~500 a）占比始终最大.排泄区硝酸盐浓度受相邻局部水流系统主导，其发育深度越大，水文遗留时间越长.若不削减输入，排泄区硝酸盐浓度经2~8 a稳定后无法下降，立即停输比逐步停输使硝酸盐提前4~6 a恢复本底值.控制地下水硝酸盐污染需优先管控相邻局部水流系统硝酸盐输入，或通过调控补给强度改变局部水流系统的发育深度，但快速削减地下水硝酸盐输入量无法显著改善管控措施对水质改善的滞后性.

Abstract

This study investigates the nitrate hydrologic residence time (NHRT) in a sandy, hierarchical groundwater flow system based on a 2D Tóth basin. Results show that with increasing groundwater recharge intensity, groundwater flow system evolves from a single regional system to nested systems, and finally to a single local flow system. The mean NHRT within the basin decreases across five distinct flow system patterns, but the proportion of longer-duration NHRT (250-500 a) remains dominant. Nitrate concentrations in discharge zones were primarily governed by adjacent local flow systems. A greater development depth of the local system prolongs NHRT in the discharge zone. Under nitrate input control scenarios, the baseline simulation (no input reduction) shows that concentrations in discharge zones stabilized within 2-8 years and then failed to decline. Compared to a gradual cessation, an immediate cessation of nitrate input accelerates recovery to background levels by only 4-6 years. The findings suggest that effective control of groundwater nitrate pollution requires either prioritizing input reduction in recharge areas of adjacent local flow systems or altering their development depth by regulating recharge intensity. Nevertheless, rapid cuts in nitrate loading do not substantially mitigate the inherent time lag of water quality response to remediation measures.

Graphical abstract

关键词

硝酸盐 / Tóth理想盆地 / 地下水流系统 / 氮遗留 / 水文遗留时间 / 污染管控措施 / 水文地质.

Key words

nitrate / Tóthian basin / groundwater flow system / legacy of nitrogen / hydrologic residence time / pollution control strategy / hydrogeology

0 引言

农业活动和大量化肥施用致使地下水硝酸盐污染日益严峻（Puckett et al.， 2011），我国28%的地下水硝酸盐含量高于世界卫生组织饮用水质量标准10 mg N/L（Gu et al.， 2013）.1961—2009年间，我国氮肥施用量（折纯量）增长了66.9倍，利用效率却不足40%（唐俊辉， 2023），未利用氮进入土壤氮库或淋滤进入含水层，引起地下水及其受纳水体硝酸盐污染.砂质含水层作为支撑农业灌溉和区域供水的关键水源地，对地表硝酸盐污染呈现出显著的环境敏感性（Spalding and Exner， 1993）.为改善地下水硝酸盐非点源污染，各类农业政策（如削减施肥量、调整耕作方式等）相继施行以控制污染物的输入，然而，土壤及地下水中遗留氮的释放会导致排泄区硝酸盐浓度的削减存在滞后性，减少补给区硝酸盐的输入短期内难以使其显著改观，甚至出现硝酸盐浓度上升的趋势（Van Meter and Basu， 2017）.因此，为达成水体硝酸盐污染削减目标，应充分考虑硝酸盐的遗留过程，以地下水系统的硝酸盐水文遗留时间（nitrate hydrologic residence time， NHRT）为依据制定管理政策.

氮肥的过量施用与其较低的利用效率导致氮素在土壤中累积（Gao et al.， 2021），在植物-土壤-微生物系统内不断运移转化，通过矿化、固定等作用产生生物地球化学遗留；经降水及农业灌溉直接淋滤、或经生物地球化学遗留转化释放后进入含水层（其中NO₃⁃N占总氮淋溶量的76.9%~88.1%），沿不同级次或长度的水流路径输移产生水文遗留.这两种氮遗留过程使得地下水系统中的历史遗留氮不断累积，并对地下水排泄区水质产生持续性的影响（胡敏鹏， 2019； 刘瀚文等， 2024）.Chen et al. （2014）在流域尺度氮均衡研究中，估算得1980—2010年储存在土壤和地下水中的总氮遗留量约为净人为总氮输入的25%.Erler et al. （2018）发现，即使研究区农业活动减少后，地下水中的遗留氮素仍在至少15 a内持续地对排泄区珊瑚产生影响.砂质沉积物较强的渗透性促使硝酸盐快速淋滤进入地下水（Chen et al.， 2018），水文遗留过程主导了硝酸盐在含水层中的停留时长.

在地下水系统中，地下水滞留时间指的是地下水从补给区运移到排泄区所经历的时间，它等于地下水年龄和预期寿命的总和（图1）.类似地，将硝酸盐从潜水面运移至排泄区或通过反硝化被完全转化为气态损失（N₂，N₂O等）而去除所经历的时间，定义为水流系统中硝酸盐的水文遗留时间（Chen et al.， 2018；胡敏鹏， 2019）.Tóth （1963）通过构建具有线性坡度、波状起伏潜水面的理想盆地模型，揭示了重力驱动下的地下水流具有局部（local， L）、中间（intermediate， M）和区域（regional， R）级次嵌套演化的自组织特征.在给定地形与渗透系数的地下水系统中，水流系统的演替主要取决于入渗补给强度的变化（Liang et al.， 2013；张萧琅和焦赳赳， 2023）.随着入渗补给强度增大，区域及中间水流系统将逐渐萎缩直至消失，转变为以局部流为主的系统模式.不同级次的地下水流往往具有不同的补给历时、相对独立的径流和排泄途径，其地下水滞留时间为数月到数万年不等（Wang et al.， 2022），其奠定了硝酸盐水文遗留时间的基本特征.同时，水流系统模式控制着含水层自补给区至排泄区的水化学分带和演化特征（张俊， 2021； Wang et al.， 2025）.危润初等（2020）对洞庭湖区地下水流系统的研究结果表明，地下水逐渐由外围的氧化环境转为湖区腹地的还原环境，硝酸盐稳定及还原带、锰（Ⅳ）还原带和铁（Ⅲ）还原带等氧化还原带依次过渡出现，直接控制着硝酸盐的生物地球化学反应，进一步影响着硝酸盐的水文遗留时间.然而，当前对多级次地下水流系统中的硝酸盐水文遗留时间特征缺乏系统的解析，相关研究鲜有报道.

本研究以小型砂质盆地为例，基于二维Tóth模型模拟计算了硝酸盐在多级次地下水流系统中的水文遗留时间，揭示水流系统模式、硝酸盐输入削减管控方式对硝酸盐水文遗留时间的影响，评估了排泄区硝酸盐恢复至本底浓度所需要的时间，为地下水硝酸盐污染防治和管理政策的制定提供科学依据.

1 研究方法

1.1　多级次地下水流系统概念模型

盆地多级次地下水流系统概念模型如图1所示.盆地地形呈波状起伏，包括多个补给区（记为RZ）和潜在排泄区（记为S）.上边界为入渗-渗流组合边界：当潜水位低于潜在排泄区高程时，地下水接受大气降水及灌溉的入渗补给，RZ处为定流量入渗边界；随着入渗补给强度增大，潜水位抬升，当潜水位高于潜在排泄区高程时，地下水从S处排泄（权董杰， 2012），此处的上边界转变为渗流边界（年庚乾等， 2020）.盆地侧缘及底部概化为隔水边界.为简化模型，将上边界的入渗补给通量设置为净入渗量，不再单独考虑蒸发过程.

降雨及灌溉水均会携带大量的溶解氧及硝酸盐经土壤入渗补给地下水，硝酸盐自地下水流系统的补给区向排泄区的运移过程中，反硝化是削减和去除氮素的主要方式.溶解态或颗粒状的有机碳是参与反硝化作用的主要电子供体（Peterson et al.， 2013），其优先被O₂氧化，阻碍电子向NO₃^-的传递，因此，当氧气几乎被完全消耗时，反硝化作用才会启动（Kolbe et al.， 2019）.电子供体对溶解氧、硝酸根等依次还原，使地下水出现相应的氧化还原分带（氧还原带、硝酸盐还原带等）.在不同级次水流系统中，水流路径和氧化还原分带的差异直接影响了硝酸盐在地下水中的遗留，且不同级次的地下水流在排泄区混合直接影响了排泄区硝酸盐的浓度.

1.2　数学模型

1.2.1　地下水流

采用余弦函数的线性组合描述图1中地面高程的变化（Jiang et al.， 2011）：

Z_{\mathrm{t}\mathrm{e}\mathrm{r}}\left(x\right)=h_{\mathrm{0}}+\left[h_{\mathrm{R}}-h_{\mathrm{R}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi x}{L}\right)\right]+\left[h_{\mathrm{L}}-h_{\mathrm{L}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{k\pi x}{L}\right)\right]，(1)

式（1）中：Z_ter为上边界地形高程（m）；h₀为盆地最低排泄点的高程（m）；L为盆地长度（m）；h_R、h_L分别为区域及局部起伏的振幅（m）；k取偶数值时，左右边界均对应山谷.

建立均质各向同性的二维剖面稳定流模型：

K\left(\frac{{\partial }^{\mathrm{2}}H}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}H}{\partial z^{\mathrm{2}}}\right)=\mathrm{0}，(2)

q_{\mathrm{1}}{|}_{{\Gamma }_{\mathrm{1}}}=I\left(p\mathrm{0}\right)\epsilon +I\left(p\ge \mathrm{0}\right)\frac{K}{l}\left(-\frac{p}{\rho g}\right)，(3)

\frac{\partial H}{\partial \mathit{n}}{\left|\right.}_{{\Gamma }_{\mathrm{0}}}=\mathrm{0}，(4)

式（2）~ （4）中：H为盆地地下水水头（m）；K为均质各向同性盆地的渗透系数（m/d）；q₁为上边界Г₁的法向通量（m/d），以进入模型域为正；I （·）为示性函数，括号内条件成立则赋值为1，反之为0；p为上边界半透水层下表面的孔隙水压力（Pa），其控制着入渗-渗流边界的转换，式3右侧第一项对应补给区入渗边界Г_RZ，第二项对应排泄区渗流边界Г_S；ε为入渗补给强度（m/d）；l为耦合长度尺度（m），代表上边界半透水层的等效厚度；K/l为半透水层对流体的传导率（1/d），表征排泄区的泄流能力；ρ为流体密度（kg/m³）；g为重力加速度（m/s²）；Г₀包括左右边界及下边界，外法线单位向量为 n .

1.2.2　反应性溶质运移

在地下水流系统补给区，硝酸盐和溶解氧浓度为定值，即定浓度边界；地下水流系统排泄区，设置为溶质自由流出边界；其余边界为零通量边界.瞬态溶质运移的数学模型如公式（5~9），对于稳态模拟，公式（5）中时间导数项为零（∂C _i /∂t = 0）.

\frac{\partial C_i}{\partial t}=\frac{\partial }{\partial j}\left(D_{jk}\frac{\partial C_i}{\partial k}\right)-\frac{\partial \left(C_i{\mathit{u}}_j\right)}{\partial j}+Q_i-\sum R_i，(5)

C_i\left(x,z\right){\left|\right.}_{t=\mathrm{0}}=C_{\mathrm{0}i}\left(x,z\right)，(6)

C_i\left(x,z,t\right){\left|\right.}_{{\Gamma }_{RZ}}=C_{\mathrm{1}i}\left(x,z,t\right)，(7)

\mathit{n}\cdot \left(-\mathit{D}\cdot \nabla C_i\right){\left|\right.}_{{\Gamma }_S}=\mathrm{0}，(8)

\mathit{n}\cdot \left(-\mathit{D}\cdot \nabla C_i+\mathit{u}{C}_i\right){\left|\right.}_{{\Gamma }_{\mathrm{0}}}=\mathrm{0}，(9)

式（5）~ （9）中：C_i 为溶质i的质量浓度（mg/L），本研究中i代表溶解氧（DO）或硝酸根（NO₃⁻）；D_jk （或 D ）为jk方向的水动力弥散系数张量（m²/d），包括机械弥散和分子扩散作用，其中j， k = x， z，在具有显著地下水流速的地下水系统（如砂质盆地等），机械弥散的贡献远大于分子扩散，因此，水动力弥散系数张量不再根据溶质i进行区分； u_j 为孔隙水流速矢量（m/d），其中j = x， z；Q_i 为溶质i的源汇项（mg/（L•d））；∑R_i 为化学反应项（mg/（L•d））；C_0i 为初始质量浓度（mg/L）；C_1i 为补给区溶质i的输入浓度（mg/L）；∇为Nabla算子； n 为边界外法向单位向量.具体地，水动力弥散系数根据如下公式进行计算：

D_{xx}={\alpha }_L\frac{u_x^{\mathrm{2}}}{\left|\mathit{u}\right|}+{\alpha }_T\frac{u_z^{\mathrm{2}}}{\left|\mathit{u}\right|}+D^{\mathrm{*}}，(10)

D_{zz}={\alpha }_T\frac{u_x^{\mathrm{2}}}{\left|\mathit{u}\right|}+{\alpha }_L\frac{u_z^{\mathrm{2}}}{\left|\mathit{u}\right|}+D^{\mathrm{*}}，(11)

D_{xz}=D_{zx}=\left({\alpha }_L-{\alpha }_T\right)\frac{u_x{u}_z}{\left|\mathit{u}\right|}，(12)

式（10）~ （12）中：α_L 和α_T 分别为土壤溶质的纵向和横向弥散度（m）；D^* 为有效分子扩散系数（m²/d）；\left|\mathit{u}\right|=\sqrt[]{u_x^{\mathrm{2}}+u_z^{\mathrm{2}}}为孔隙水流实际流速的大小（m/d）.

砂质沉积物中固氮作用和有机氮矿化作用较弱，反硝化作用是去除地下水中硝态氮的主要过程，外源有机质的持续入渗可为氧气和硝酸盐的还原过程提供丰富的电子供体（刘晨等， 2025），使其反应速率主要受氧气及硝酸盐浓度控制，符合一级反应动力学特征（Green et al.， 2016）：

R_{\mathrm{D}\mathrm{O}}=k_{\mathrm{D}\mathrm{O}}{C}_{\mathrm{D}\mathrm{O}}，(13)

R_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}=I\left(C_{\mathrm{D}\mathrm{O}}\mathrm{2}\mathrm{m}\mathrm{g}/L\right)k_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}{C}_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}，(14)

式（13）~ （14）中：k_DO、k{}_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}分别为氧气和硝酸盐还原速率常数（a^-1）；C_DO、C{}_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}为t时刻两种溶质的浓度（mg/L）；I （·）为示性函数.

1.2.3　地下水滞留时间

盆地地下水滞留时间是地下水年龄与预期寿命的总和（图1），体现了地下水资源的更新能力.Goode（1996）引入“年龄质量”的概念，将地下水年龄类比为地下水中运动的溶质，采用对流弥散方程对其进行描述，得到稳定状态下地下水年龄控制方程（式15）.其反映了地下水的理论年龄，为分析盆地中地下水年龄的空间分布规律提供了便利（蒋小伟 等， 2012）.

\nabla \cdot \left(\theta \mathit{D}\cdot \nabla \tau \right)-\nabla \cdot \left(\mathit{u}\theta \tau \right)+\theta =\mathrm{0}，(15)

式（15）中：τ为地下水平均年龄（a）；θ为孔隙度.

地下水年龄为地下水质点从地下水面运动到当前位置所需时间，而地下水预期寿命等于地下水质点从当前位置运移到最终排泄出盆地所要经历的时间（王俊智， 2015），其计算方法与式（15）形式相似，只是将对流项取相反数.

\nabla \cdot \left(\theta \mathit{D}\cdot \nabla {\tau }_{\mathrm{e}\mathrm{x}\mathrm{p}}\right)-\nabla \cdot \left(\mathit{u}\theta {\tau }_{\mathrm{e}\mathrm{x}\mathrm{p}}\right)+\theta =\mathrm{0}，(16)

式（16）中：τ_exp为地下水预期寿命（a）.

在水流模型的隔水边界上，年龄及预期寿命的质量通量 uτ⁃ D ·\nablaτ均为0.对于地下水年龄，上边界的补给区为已知浓度边界，即τ为0，排泄区为对流边界，即弥散通量 D ·\nablaτ为0；对于地下水预期寿命，上边界的排泄区为已知浓度边界，即τ_exp为0，补给区为对流边界，弥散通量 D ·\nablaτ_exp为0.

1.3　硝酸盐水文遗留时间计算方法

本文重点研究砂质含水层中硝酸盐的水文遗留过程，暂不考虑土壤中氮的固定等生物地球化学遗留过程.硝酸盐的水文遗留时间指其在水流系统中停留时间，即硝酸盐从地下水面运移至排泄区势汇或通过反硝化气态损失（N₂，N₂O等）被完全去除所经历时间.结合水流系统理论，本文自主构建以下计算方法：对于能完整穿越地下水系统的硝酸盐质点，其水文遗留时间与地下水滞留时间相同；对于不能完整穿越地下水系统、中途即被完全反硝化的硝酸盐质点，其水文遗留时间等于被完全反硝化所需时间，即溶解氧及硝酸盐还原的总用时.以硝酸根浓度C{}_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}_{， threshold}=0.02 mg/L为临界阈值，当硝酸盐浓度低于该值时，视为已达到完全去除的下限标准（图1）.

t_{\mathrm{N}\mathrm{H}\mathrm{R}\mathrm{T}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(t_{\mathrm{D}\mathrm{O}}+t_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}},\tau +{\tau }_{\mathrm{e}\mathrm{x}\mathrm{p}}\right)，(17)

t_{\mathrm{D}\mathrm{O}}=\frac{\mathrm{1}}{k_{\mathrm{D}\mathrm{O}}}\mathrm{l}\mathrm{n}\frac{C_{\mathrm{D}\mathrm{O},\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}}}{C_{\mathrm{D}\mathrm{O},\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}}}，(18)

t_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}=\frac{\mathrm{1}}{k_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}}\mathrm{l}\mathrm{n}\frac{C_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}},\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}}}{C_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}},\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}}}，(19)

式（17）~ （19）中：t_NHRT为硝酸盐水文遗留时间（a）；t_DO为溶解氧消耗的用时（a），即反硝化滞后时间；t{}_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}为硝酸盐还原的反应用时（a）；τ+τ_exp为地下水滞留时间（a）；C_{DO， watertable}、C{}_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}_{， watertable}分别为溶解氧、硝酸根的输入浓度（mg/L）；C_{DO， threshold}、C{}_{\mathrm{N}{\mathrm{O}}_{\mathrm{3}}}_{， threshold}分别为反硝化启动时的DO阈值浓度及硝酸根完全去除的阈值浓度（mg/L）.

由于水流路径的差异，同一时刻进入盆地的硝酸盐质点在盆地中的水文遗留时间并不相同.基于欧拉法模拟结果，在所有满足硝酸盐检出限要求的网格单元中，获取大量NHRT样本数据.通过构建概率密度函数，定量表征水流系统硝酸盐水文遗留时间的分布特征.NHRT分布集总了地下水流场的时空信息，反映了水流系统中地下水遗留硝酸盐的汇集情况.

而在水流系统排泄区，根据硝酸盐浓度随时间的变化，将排泄区硝酸盐的水文遗留时间t'_NHRT定义为干预（减少）地下水硝酸盐的输入时刻t_c至排泄区恢复到硝酸盐本底值所需的时间（式20），该值反映了控制硝酸盐输入的管控实施下，地下水遗留氮库中的硝酸盐完全削减所需要的时间（Gao et al.， 2021）.

t_{\mathrm{N}\mathrm{H}\mathrm{R}\mathrm{T}}^{\mathrm{\text{'}}}=t_{\mathrm{b}\mathrm{g}}-t_{\mathrm{c}}，(20)

式（20）中：t_bg为管控措施施行后，排泄区硝酸盐浓度首次下降到环境本底浓度后的时刻.

2 模型设置

2.1　模拟情景及参数

参照Udden⁃Wentworth粒级分类标准及相关经验数据（Wentworth， 1922； Allard and World Health Organization， 1992； Domenico and Schwartz， 1997），砂质沉积物粒径通常大于64 μm，其渗透系数K与孔隙度θ的经验范围分别为1.72×10^-2~518.4 m/d、0.3~0.4.本研究模拟盆地宽度L为3 600 m，长宽比约为10∶1，含水介质为均质各向同性的砂质沉积物，K取值0.25 m/d，θ取值0.40（表1），主要对应细砂介质.溶质的弥散度具有尺度效应，纵向弥散度α_L约取盆地宽度L的1/100，横向弥散度α_T取纵向弥散度α_L的1/10（Gelhar et al.， 1992），分子扩散系数D^*取1.16×10^-9 m²/s.假定包气带的水分中溶解氧与空气达到饱和，以该饱和平衡的溶解氧浓度（C_{DO， watertable} =6.4 mg/L）作为模型上边界的溶解氧输入浓度（孙锐， 2018）.地下水溶解氧运移过程中与电子供体反应而衰减，当溶解氧浓度低于C_{DO， threshold}=2 mg/L时（Korom， 1992；蒋小伟 等， 2012），电子供体在微生物介导下参与反硝化作用.McMahon et al. （2008）估计了美国Nebraska砂质潜水含水层中氧气还原和反硝化反应的一级反应动力学常数，分别介于0.02~0.30 a^-1和未检出（n.d.）~0.05 a^-1范围内.为模拟砂质潜水含水层中DO和NO₃^-的一级动力学衰减过程，取k_DO=0.04 a^-1，k\mathrm{N}{\mathrm{O}}_{\mathrm{3}}=0.02 a^-1，且不随时空变化.

通过改变上边界入渗补给强度（ε），即调整入渗强度比R_ic （R_ic = ε/K），在稳态条件下模拟对应的5种水流系统模式，如表2所示.

为探究水流系统模式、硝酸盐输入削减管控方式对硝酸盐水文遗留时间的影响，设置3组工况.工况1（Case⁃I）对比稳态条件下5种水流系统模式（R、RL、RML、LM、L）中NHRT分布特征；基于Case⁃I中5种稳态水流系统模式，工况2、3（Case⁃II和Case⁃III）分别探究硝酸盐输入时空削减管控措施（硝酸盐输入位置、输入强度）对不同水流系统排泄区NHRT的影响，通过预测排泄区硝酸盐的浓度变化，评估水体恢复至环境本底值所需时间.

3组工况的水流模型设置均相同，区别在于溶质运移模型中硝酸盐的时空输入方式（表3）.其中，Case⁃II设置3组空间管控情景，在不同补给区（RZ1、RZ2、RZ3）输入硝酸盐污染物，待硝酸盐运移稳定后（浓度年际变化低于0.2%）停止输入，对比不同级次水流控制下的排泄区NHRT特征；Case⁃III中，硝酸盐的输入浓度参照我国1961—2020年间单位面积农业区的氮肥淋失量（国际粮农组织，https：//www.fao.org/faostat）的变化趋势，以2020年为现状年，对比3种地下水硝酸盐输入强度管控效果（图2）：情景（ⅰ）“保持现状”，维持现有的土地利用及管理措施，假设自2020年起硝酸盐输入保持在当前值；情景（ⅱ）“逐步停止”，参照全国氮肥淋失量的变化，假设淋失进入地下水的硝酸盐以每年20%的速率减少，直至降至0 mg/L；情景（ⅲ）“立即停止”，假设立即停止硝酸盐输入地下水.

2.2　模型求解

利用COMSOL Multiphysics 6.3求解地下水中硝酸盐的迁移转化过程及水文遗留时间.运用Delaunay三角网格剖分将几何模型进行网格离散化，并加密盆地边界附近的剖分单元，网格平均单元质量0.902，共剖分有效网格10⁴余个.Case⁃II、Case⁃III模拟期均为250 a，为提高计算效率，采用非均匀时间步长.硝酸盐输入变化较大的时间段（Case⁃II第81~89年、Case-III中2021—2029年）取时间步长为1 a，Case⁃II其余时刻时间步长为5 a，Case⁃III在2100年后取时间步长为2 a，其余时刻时间步长均为5 a.

3 结果与讨论

3.1　Case⁃I：不同水流系统模式下的硝酸盐水文遗留时间分布特征

受地形和地质条件控制，随入渗补给强度的增大，盆地潜水位抬升，局部水流系统发育深度及数量增加，而区域水流系统、中间水流系统逐渐消失，依次呈现5种水流系统模式（图3Aa~3Ae）：单一区域水流系统（R）→局部-区域二级嵌套水流系统（2L+R）→局部-中间-区域三级嵌套水流系统（6L+M+R）→局部-中间二级嵌套水流系统（6L+M）→单一局部水流系统（6L）.

溶解氧及硝酸盐的模拟浓度稳定后，呈现出硝酸盐浓度沿着地下水渗流路径不断降低、硝酸盐还原带的宽度随入渗补给强度的增大而减小的特点（图3B）.单一区域水流系统（R，图3Ba）中，地下水氧还原带仅分布在最强势汇处，最大深度不足15 m，含水层中较大的缺氧（低于2 mg/L）面积有利于地下水体中反硝化作用的发生；局部-中间-区域三级水流系统（6L+M+R，图3Bc）中，地下水硝酸盐还原范围宽为230~270 m；单一局部水流系统（6L，图3Be）中，地下水硝酸盐还原范围仅宽为200~250 m.硝酸盐还原带的范围的不断缩减，阻碍了反硝化作用对地下水硝酸盐的削减，导致更多的硝酸盐在地下水系统中遗留.

盆地地下水滞留时间分布与地下水流系统的发育密切相关.受径流路径与地下水流速的影响，区域水流系统中地下水滞留时间最长，中间水流系统次之，局部水流系统的滞留时间最短.在不同的水流模式中，随入渗补给强度的增大，地下水平均滞留时间逐渐减小（表4）.

相较地下水滞留时间，NHRT整体较短，且随局部水流系统发育深度的增大，其分布逐渐离散，较长NHRT（250~500 a）的占比始终最大（表4、图4）.单一区域水流系统（R）中，遗留的硝酸盐所占空间范围仅为地下水系统的13%（图4Aa），但平均NHRT最长（455 a）；局部-区域二级嵌套水流系统（2L+R）中，硝酸盐主要分布在局部水流系统L1、L2中（图4Ab），平均NHRT为452 a；局部-中间-区域三级嵌套水流系统（6L+M+R）中，硝酸盐在局部和中间水流系统的全部区域均有分布（图4Ac），平均NHRT为421 a；在前述三种水流系统模式中，硝酸盐反硝化削减量超过80%（表4）.相较地下水滞留时间，NHRT更短且分布更集中（图4Ba~Bc），即包含硝酸盐的“老水”所占比例较低，表明较长的运移路径及缓慢的流速有利于反硝化作用的充分进行.

在局部-中间二级嵌套水流系统（6L+M）、单一局部水流系统（6L）中，反硝化对硝酸盐的削减量分别为76%、68%，尚未完全反硝化的硝酸盐几乎遍布整个地下水系统（图4Ad、4Ae），NHRT与地下水滞留时间的分布特征愈加接近（图4Bd、4Be）.局部水流系统硝酸盐的还原能力较弱，不存在完全反硝化的区域，而随着入渗补给强度的增大，中间水流系统、区域水流系统逐渐衰退，反硝化的能力显著下降，硝酸盐迁移及污染范围不断增大.较高的地下水补给通量一方面导致还原环境的范围减小，抑制反硝化作用的发生，另一方面，地下水流速的加快及径流路径的缩短，减少了地下水携带的NO₃^-的水力停留时间，更多的硝酸盐不经过反硝化作用直接到达排泄区（Liu et al.， 2017；王佳琪 等， 2019）.

3.2　Case⁃II：不同输入位置对排泄区硝酸盐水文遗留时间的影响

Case⁃II模拟结果显示，在局部-中间-区域水流系统中，补给区RZ1处输入的硝酸盐可分别通过局部水流系统L1、L2输移至排泄区S1、排泄区S2（图3Ac），对应的硝酸盐稳定浓度分别为29 mg/L、10 mg/L；硝酸盐输入停止后，S1、S2对应的NHRT分别为135 a、50 a（图5Ac）；而在补给区RZ3输入时，S1、S2均无响应（≤ 1 mg/L）（图5Cc）.表明本例中局部水流系统对硝酸盐的削减能力有限，硝酸盐不能完全反应，而区域水流系统、中间水流系统因较慢的运移速度及较强的还原环境，反硝化作用较为充分，与Case⁃I呈现的不同级次水流系统中硝酸盐的削减规律一致.相较局部水流系统L2，L1的发育深度更大，在含水层深部，未通过反硝化作用去除的硝酸盐，对排泄区贡献了大量硝酸盐负荷，同时延长了排泄区硝酸盐的水文遗留时间.当水流系统模式不变时，发育深度越大的局部水流系统对应的排泄区硝酸盐稳定浓度越大、停止硝酸盐输入后NHRT越长.此外，浅部极短水文路径中，相对较快的水循环速度及微弱的反硝化作用使得单一区域水流系统中的排泄区S1也表现出较高的硝酸盐浓度及较短的NHRT（图5Aa），局部-区域水流系统中的L2、各水流模式中的L6所对应的排泄区（图3、图5Bb、图5Cc~5Ce）也表现出同样的情况.

当仅RZ1或RZ2输入硝酸盐时，随着入渗补给强度的增大（图5Ab~5Ae、5Bc~5Be），水流系统模式由局部-区域嵌套系统向单一局部水流系统转化，在发育深度更大的局部水流系统中（如L1、L3），排泄区硝酸盐的稳定浓度增长趋缓，硝酸盐浓度拖尾减弱，NHRT逐渐缩短；而在发育深度较小的局部水流系统中（如L2、L4），排泄区硝酸盐稳定浓度的增幅较大，NHRT迅速增加（图5Bc~5Be）.当水流系统模式发生改变时，不同发育深度的局部水流系统中，排泄区NHRT呈现出两种演变规律：在发育深度较大的局部水流系统中，地下水运移深度的进一步增加使得反硝化反应更为充分，反硝化作用缩短了排泄区硝酸盐的水文遗留时间；而在发育深度较小的局部水流系统中，反硝化作用的增强有限，但更长的运移路径导致硝酸盐遗留时间增加.

3.3　Case⁃III：不同输入强度对排泄区硝酸盐水文遗留时间的影响

Case⁃III中，各排泄区硝酸盐的水文遗留时间演变规律与Case⁃II相似，即随着局部水流系统发育程度增大，NHRT逐渐趋同并整体减小.但Case⁃III在各补给区同时输入硝酸盐，使得排泄区受到多个相邻局部水流系统流出通量的影响，从而缩小了各排泄区硝酸盐的浓度差异（图6）.“保持现状”情景中（图6Aa~6Ae），不同水流系统模式的排泄区硝酸盐浓度逐渐升高，并后均在管控实施后的2~8 a内趋于稳定（浓度变化0.1%）.单一局部水流系统（R）中，硝酸盐的输出浓度约在管控措施实行8 a后（2028年）达到稳定（28 mg/L），反硝化作用削减约72%的硝酸盐；局部-中间-区域三级水流系统中，硝酸盐的输出浓度约在管控措施实行4~8 a后（2024—2028年）达到稳定（31~75 mg/L）；单一局部水流系统中，各排泄区稳定浓度（54 mg/L）及达稳定所需时间（2~6 a）差异不大，反硝化作用仅削减46%的硝酸盐.表明“保持现状”情景中，局部水流系统发育深度的增大不利于削减硝酸盐浓度，排泄区硝酸盐污染程度加剧.

在两种减少硝酸盐输入的管控措施实施后，水流系统排泄区需94~142 a恢复至环境本底值.相较“逐步停止”情景，“立即停止”情景中排泄区硝酸盐浓度在管控实施初期（约前10 a）下降速度更快，但NHRT并未显著减少（图6Ba~6Be、图6Ca~6Ce）.在各水流系统模式中，硝酸盐输入“立即停止”时排泄区硝酸盐提前恢复至环境本底值的时间优势仅4~6 a.

4 结论

本文以小型砂质潜水盆地为例，模拟计算了不同水流系统模式下的硝酸盐水文遗留时间，评估了不同硝酸盐输入管控措施下排泄区地下水硝酸盐浓度恢复至本底值所需的时间，得到以下结论.

（1）随着入渗补给强度增大，依次呈现单一区域水流系统（R）、局部-中间-区域三级嵌套水流系统（6L+M+R）、单一局部水流系统（6L）等5种水流系统模式，反硝化对硝酸盐的削减比例随之逐渐降低，硝酸盐迁移及污染范围从整个模型区域的13%扩大至98%，NHRT均值从455 a逐渐降至336 a，分布也愈加集中，逐渐趋近于地下水滞留时间分布，而受较长运移路径的影响，较长水文遗留时间（250~500 a）的占比始终最大.

（2）地下水排泄区的硝酸盐浓度主要由相邻局部水流系统控制，受区域水流系统、中间水流系统的影响较小（≤1 mg/L）.当水流系统模式不变时，局部水流系统发育深度越大，排泄区硝酸盐的浓度峰值越大、水文遗留时间越长.当水流模式由局部-区域嵌套系统向单一局部水流系统转化时，排泄区NHRT呈现不同的演变规律：在发育深度较大的局部水流系统中，增强的反硝化作用主导着NHRT的缩短现象；而发育深度较小的局部水流系统中，延长的地下水渗流路径及地下水滞留时间，反而导致NHRT进一步增长.

（3）若不减少硝酸盐污染输入，各水流系统模式中排泄区的硝酸盐浓度约在管控施行后的2~8 a内达到稳定（浓度变化0.1%）后不再下降，反硝化作用对硝酸盐的削减比例为25%~83%.相较逐步停止硝酸盐输入，立即停止硝酸盐输入仅能缩短排泄区硝酸盐的水文遗留时间4~6 a.

为减少地下水硝酸盐污染，应首先控制相邻局部水流系统补给区的硝酸盐输入，或通过调控补给强度改变局部水流系统的发育深度.若不减少人为硝酸盐输入而仅依靠地下水系统中的生物地球化学作用削减硝酸盐负荷，局部水流系统的排泄区仍将产生硝酸盐污染问题.尽快削减地下水硝酸盐输入量在管控实施初期可迅速降低排泄区硝酸盐浓度，但不能显著改善污染管控措施对水质改善的滞后性.实际工程中，可结合水文地质调查与遥感分析手段，识别地形驱动下的潜在势源与势汇；利用环境同位素和测年等手段识别滞留时间较短的局部水流路径，圈定关键补给区，进而通过精准控肥或调整土地利用方式，对多级次水流系统中的硝酸盐污染实现源头管控.

本文基于二维Tóth模型探究了硝酸盐在多级次地下水流系统中的水文遗留时间特征，但在刻画自然界复杂地质条件方面仍存在局限.未来应进一步考虑介质非均质性、气候瞬态变化及地表水-地下水交互过程等影响，以期为区域尺度地下水环境修复与动态管理提供更精准的预测依据.

参考文献

[1] Allard, M., World Health Organization, 1992. GEMS/WATER Operational Guide. Report GEMS/W.92.1. National Water Research Institute, Burlington.

[2] Chen, D. J., Hu, M. P., Dahlgren, R. A., 2014. A Dynamic Watershed Model for Determining the Effects of Transient Storage on Nitrogen Export to Rivers. Water Resources Research, 50(10): 7714-7730. https://doi.org/10.1002/2014wr015852

[3] Chen, D. J., Shen, H., Hu, M. P., et al., 2018. Legacy Nutrient Dynamics at the Watershed Scale: Principles, Modeling, and Implications. Advances in Agronomy. Elsevier,Amsterdam,237-313. https://doi.org/10.1016/bs.agron.2018.01.005

[4] Domenico, P. A., Schwartz, F. W., 1997. Physical and Chemical Hydrogeology. John Wiley Sons, New York.

[5] Erler, D. V., Shepherd, B. O., Linsley, B. K., et al., 2018. Coral Skeletons Record Increasing Agriculture-Related Groundwater Nitrogen Inputs to a South Pacific Reef over the Past Century. Geophysical Research Letters, 45(16): 8370-8378. https://doi.org/10.1029/2018gl078656

[6] Gao, J. B., Wang, S. M., Li, Z. Q., et al., 2021. High Nitrate Accumulation in the Vadose Zone after Land⁃Use Change from Croplands to Orchards. Environmental Science Technology, 55(9): 5782-5790. https://doi.org/10.1021/acs.est.0c06730

[7] Gelhar, L. W., Welty, C., Rehfeldt, K. R., 1992. A Critical Review of Data on Field⁃Scale Dispersion in Aquifers. Water Resources Research, 28(7): 1955-1974. https://doi.org/10.1029/92wr00607

[8] Goode, D. J., 1996. Direct Simulation of Groundwater Age. Water Resources Research, 32(2): 289-296. https://doi.org/10.1029/95wr03401

[9] Green, C. T., Jurgens, B. C., Zhang, Y., et al., 2016. Regional Oxygen Reduction and Denitrification Rates in Groundwater from Multi⁃Model Residence Time Distributions, San Joaquin Valley, USA. Journal of Hydrology, 543: 155-166. https://doi.org/10.1016/j.jhydrol.2016.05.018

[10] Gu, B. J., Ge, Y., Chang, S. X., et al., 2013. Nitrate in Groundwater of China: Sources and Driving Forces. Global Environmental Change, 23(5): 1112-1121. https://doi.org/10.1016/j.gloenvcha.2013.05.004

[11] Hu, M. P., 2019. Quantitative Study on Lag Effect of Watershed Non⁃Point Source Nitrogen Pollution (Disserattion). Zhejiang University, Hangzhou (in Chinese with English abstract).

[12] Jiang, X. W., Wan, L., Wang, X. S., et al., 2012. Distribution of Groundwater Age in Drainage Basins. Hydrogeology Engineering Geology, 39(4): 1-6 (in Chinese with English abstract).

[13] Jiang, X. W., Wang, X. S., Wan, L., et al., 2011. An Analytical Study on Stagnation Points in Nested Flow Systems in Basins with Depth⁃Decaying Hydraulic Conductivity. Water Resources Research, 47: 2010WR009346. https://doi.org/10.1029/2010wr009346

[14] Kolbe, T., de Dreuzy, J. R., Abbott, B. W., et al., 2019. Stratification of Reactivity Determines Nitrate Removal in Groundwater. Proceedings of the National Academy of Sciences of the United States of America, 116(7): 2494-2499. https://doi.org/10.1073/pnas.1816892116

[15] Korom, S. F., 1992. Natural Denitrification in the Saturated Zone: A Review. Water Resources Research, 28(6): 1657-1668. https://doi.org/10.1029/92wr00252

[16] Liang, X., Quan, D. J., Jin, M. G., et al., 2013. Numerical Simulation of Groundwater Flow Patterns Using Flux as Upper Boundary. Hydrological Processes, 27(24): 3475-3483. https://doi.org/10.1002/hyp.9477

[17] Liu, C., Gong, X. L., Liang, Y., et al., 2025. Characteristics of Seasonal Changes in Organic Matter of Groundwater in Binhai, Jiangsu Province and Its Impact on Nitrogen Transport and Transformation. Earth Science, 50(6): 2400-2415 (in Chinese with English abstract).

[18] Liu, H. W., Luan, H. A., Zhang, Y. T., et al., 2024. Research Progress on the Occurrence Characteristics of Nitrogen and Phosphorus Leaching in Dryland Farmland. Chinese Journal of Eco⁃Agriculture, 32(9): 1520-1533 (in Chinese with English abstract)

[19] Liu, Y. Y., Liu, C. X., Nelson, W. C., et al., 2017. Effect of Water Chemistry and Hydrodynamics on Nitrogen Transformation Activity and Microbial Community Functional Potential in Hyporheic Zone Sediment Columns. Environmental Science Technology, 51(9): 4877-4886. https://doi.org/10.1021/acs.est.6b05018

[20] McMahon, P. B., Böhlke, J. K., Kauffman, L. J., et al., 2008. Source and Transport Controls on the Movement of Nitrate to Public Supply Wells in Selected Principal Aquifers of the United States. Water Resources Research, 44(4): 2007WR006252. https://doi.org/10.1029/2007wr006252.

[21] Nian, G. Q., Chen, Z. H., Zhang, L. F., et al., 2020. Treatment of Two Boundary Conditions for Rainfall Infiltration in Slope and Its Application. Rock and Soil Mechanics, 41(12): 4105-4115 (in Chinese with English abstract).

[22] Peterson, M. E., Curtin, D., Thomas, S., et al., 2013. Denitrification in Vadose Zone Material Amended with Dissolved Organic Matter from Topsoil and Subsoil. Soil Biology and Biochemistry, 61: 96-104. https://doi.org/10.1016/j.soilbio.2013.02.010.

[23] Puckett, L. J., Tesoriero, A. J., Dubrovsky, N. M., 2011. Nitrogen Contamination of Surficial Aquifers—A Growing Legacy. Environmental Science Technology, 45(3): 839-844. https://doi.org/10.1021/es1038358.

[24] Quan, D. J., 2012. Analysis of Two⁃Dimension Groundwater Flow Patterns and Transformation Rules (Dissertation). China University of Geosciences (Wuhan)，Wuhan(in Chinese with English abstract).

[25] Spalding, R. F., Exner, M. E., 1993. Occurrence of Nitrate in Groundwater: A Review. Journal of Environmental Quality, 22(3): 392-402. https://doi.org/10.2134/jeq1993.00472425002200030002x.

[26] Sun, R., 2018. Numerical Simulation of Coastal Nutrient Input Processes Using Transport⁃Reaction Numerical Model (Dissertation). China University of Geosciences (Beijing)，Beijing (in Chinese with English abstract).

[27] Tang, J. H., 2023. Effects of Nitrogen Reduction Combined with Urease/Nitrification Inhibitors on Nitrogen Loss and Nitrogen Use Efficiency in Purple Soil⁃Wheat System (Dissertation). Sichuan Agricultural University，Chengdu (in Chinese with English abstract).

[28] Tóth, J., 1963. A Theoretical Analysis of Groundwater Flow in Small Drainage Basins. Journal of Geophysical Research, 68(16): 4795-4812. https://doi.org/10.1029/jz068i016p04795.

[29] Van Meter, K. J., Basu, N. B., 2017. Time Lags in Watershed⁃Scale Nutrient Transport: An Exploration of Dominant Controls. Environmental Research Letters, 12(8): 084017. https://doi.org/10.1088/1748⁃9326/aa7bf4

[30] Wang, J. L., Jin, M. G., Jia, B. J., et al., 2022. Numerical Investigation of Residence Time Distribution for the Characterization of Groundwater Flow System in Three Dimensions. Journal of Earth Science, 33(6): 1583-1600. https://doi.org/10.1007/s12583⁃022⁃1623⁃3

[31] Wang, J. Q., Ma, R., Sun, Z. Y., 2019. Reactive Transport and Model of Nitrogen Pollutants in the Surface Water⁃Ground Water Interaction Zones: A Review. Geological Science and Technology Information, 38(4): 270-280 (in Chinese with English abstract).

[32] Wang, J. Z., 2015. A Methodological Study on the Identification of Hierarchically Nested Groundwater Flow Systems in Drainage Basins (Dissertation). China University of Geosciences (Beijing), Beijing(in Chinese with English abstract).

[33] Wang, Y. A., Zhao, R. C., Ding, L., et al., 2025. Karst Trough Control of Solute Transport Processes at Two Karst Groundwater Flow Systems, Western Hubei, Central China. Journal of Earth Science, 36(4): 1731-1741. https://doi.org/10.1007/s12583⁃022⁃1665⁃6

[34] Wei, R. C., Tang, S. M., Wu, C. S., et al., 2020. Redox Zoning of Shallow Groundwater in Dongting Lake Region. China Environmental Science, 40(4): 1715-1722 (in Chinese with English abstract)

[35] Wentworth, C. K., 1922. A Scale of Grade and Class Terms for Clastic Sediments. The Journal of Geology, 30(5): 377-392. https://doi.org/10.1086/622910.

[36] Zhang, J., 2021. Hierarchically Nested Groundwater Flow System and Hydrological and Ecological Effects in the Northern Ordos Basin (Dissertation). Chang’an University，Xi’an(in Chinese with English abstract).

[37] Zhang, X. L., Jiao, J. J., 2023. Numerical Modelling Study on Non⁃Steady⁃State Groundwater Flow Systems in Response to Changing Rainfall. Bulletin of Geological Science and Technology, 42(4):154-161,169 (in Chinese with English abstract).

基金资助

国家自然科学基金青年项目(42407075)

国家资助博士后研究人员计划项目(GZC20241601)

中国博士后科学基金项目(2025M773158)

自然资源部2024年度部省合作项目(2024ZRBSHZ021)