1 引言(Introduction)

基流一般定义为“来源于地下水及其他滞后水源的径流”(Hall, 1968), 是河川径流在枯水季节的基本流量, 也是河川径流在丰水季节的重要补给来源, 对流域水资源利用和生态环境保护具有非常重要的作用(He et al., 2016；于艺鹏等, 2016).与地表直接径流相比, 由于受到地下水库的调蓄, 基流具有稳定性和滞后性.从污染物源解析的角度分析, 地表直接径流和基流都可以作为流域养分进入河流水体的水文通道(薛禹群等, 2009；He et al., 2016).相应地, 河流养分负荷的来源也可分为地表直接径流源和基流源.考虑到地下水水文循环的滞后性, 许多流域, 尤其是大流域, 地下水中的养分可能是数年、甚至数十年前累积的历史产物(Meter et al., 2015).因此, 地下水以基流形式输出的养分往往具有滞后性, 无法与地表直接径流养分“同步同源”入河.换句话说, 研究时期内河流养分负荷可能有相当大的比例来源于过去数年甚至几十年的氮输入.

鉴于基流与地表直接径流养分输出在时间和途径上的巨大差异, 量化基流养分输出对于完善水源地水污染控制策略具有重要的科学意义.第一, 只有明确基流的贡献才能够将历史的污染责任分离出来, 避免将现今地表水的污染都归咎于现时污染源, 高估当前污染物排放水平, 导致严重误导污染源解析结果；第二, 由于地下水污染改善与流域污染治理措施存在长达数年以上的响应时差(Sanford et al., 2013; Vervloet et al., 2018), 在基流养分输出较高的流域, 水质的进一步提升就会成为一个较长期的问题, 需要正确引导政府和公众的预期；第三, 在水污染治理过程中, 需要充分考虑基流和地表直接径流养分入河途径上的不同, 根据二者的贡献大小, 搭配使用不同的源汇控制策略.例如, 对于地表直接径流污染贡献率较高的流域, 可实施“养分源头控制”、“河道过程拦截”、“陆地过程拦截”等措施.反之, 则应该主要采用“河道过程拦截”为主.

随着非点源污染研究的逐步深入, 围绕基流养分负荷定量问题, 学者们通过构建统计模型等手段, 已经进行了一些探索, 其中, 最为常用的是由Schilling等(2004)提出的方法.该方法以基流流量分割为基础, 将总径流中基流流量大于90%时的养分浓度近似地视为基流中养分的浓度；然后, 分别建立河流养分负荷量及基流(不含地表直接径流)养分负荷量的Load Estimator(LOADEST)/统计模型；最后, 计算出河流逐日养分负荷量和基流养分负荷量.Schilling等(2004)基于美国Raccoon河长达28年的流量和硝态氮监测数据, 利用该方法计算出基流的硝态氮输出通量占Raccoon河年平均硝态氮输出通量的66%左右.同样地, 李文超等(2018)以位于云南省的高原湖泊典型流域-凤羽河流域为例, 通过该方法揭示了基流是高原湖泊流域总氮(TN)输出的主要途径, 基流输出TN负荷占流域TN输出负荷的71.1%.除了上述方法之外, Miller等(2016)利用基流分割方法与连续高频率的硝态氮监测数据(1次·15 min^-1)建立了硝态氮浓度与基流指数的回归方程, 估算了基流与地表直接径流的硝态氮浓度, 进而计算出美国Chesapeake湾流域基流对地表水硝态氮入河负荷量的贡献率为58%~73%.Kim等(2010)在基流分割的基础上, 结合地下水硝态氮监测数据, 估算了韩国Daejeon市基流的硝态氮输出通量, 结果表明约60%的地表水硝态氮来自于流域基流的输出.然而, Schilling等(2004)和Miller等(2016)的方法对数据依赖程度高, 要求有长时间序列或高频率监测的水文水质数据, 否则无法挑选出足够的基流时期样本进行后续建模, 难以在监测体系不完善的地区推广和使用.Kim等(2010)仅以一个地下水监测井的数据代表基流浓度, 未考虑地下水中养分的时空异质性及其入河过程中的迁移转化, 研究结果可能存在较大误差.

饮用水水库上游河流大多属于山溪性河流, 其水量受基流补给强度较大, 故基流有可能是流域氮素输出的主要方式.然而, 目前关于水源地基流养分输出贡献方面的研究鲜有报道.更值得注意的是, 我国饮用水水源地, 尤其是农村分散式水源地, 监测体系还有待于完善, 基础数据较为缺乏.基于以上考虑, 本研究尝试建立一种建模简单、数据依赖程度低的水源地基流负荷分割方法.选择位于浙江省珊溪水源地的玉泉溪流域为研究区域, 基于2010—2013年逐月的TN监测数据和逐日的流量数据, 利用LOADEST模型模拟河流段末逐日TN负荷量之后, 通过改进数字滤波算法中滤波参数的计算方法进而使其能够应用于基流负荷分割, 实现了水源地基流对河流TN负荷量贡献的定量识别, 以期为水源地水质进一步改善提供科学依据.

2 研究区域(Study area)

珊溪水源地位于浙江省温州市飞云江上游, 其地理位置介于东经119°47′~120°15′和北纬27°36′~ 27°50′之间, 流域面积为2302 km²(图 1).该水源地由珊溪水库、赵山渡水库和配套引水渠系3部分组成, 控制着温州地区将近70%的水资源, 供水范围涉及到当地500多万人口, 占浙江省水库供水人口的20%以上.流域属江南丘陵山地红壤区, 土壤主要是红壤, 其余为黄壤、紫色土、粗骨土、潮土和水稻土(梅琨等, 2016).水源地区域多年(1950—2013年)平均降水量为1880 mm, 平均降水天数为149 d, 降雨量年内分配不均匀, 75%的降水集中在4—9月；平均气温、总日照时数和相对湿度分别为19.6 ℃, 1887时和83%.珊溪水源地入库支流较多, 其中, 玉泉溪是距离水源地集水区出口(赵山渡水库取水口)最近的主要支流, 水库供水水质达标受其养分输入影响较大.因此, 选择玉泉溪流域为研究对象(图 1).玉泉溪流域面积约为264 km², 流域内土地利用类型主要分为林地、耕地、建设用地和未利用地, 分别占流域总面积的64.9%、32.1%、1.8%和1.2%(图 1).自2007年以来, 随着库区环境整治工作的推进, 玉泉溪流域内工业企业基本上都已外迁, 但由于河流水质受到流域内非点源养分输出的影响, 水质状况仍不乐观.根据《珊溪(赵山渡)库区水环境综合整治和生态保护规划(2011—2020年)》和温州市珊溪水利枢纽管理局水质监测数据显示, 2010—2013年玉泉溪入库断面未达到稳定符合地表水Ⅱ类水质的要求, TN是主要超标因子.

3 数据与方法(Data and methods) 3.1 数据

2010—2013年间的玉泉溪段末水质监测数据和珊溪水库流量数据由温州市珊溪水利枢纽管理局提供, 水质监测频率为1次·mon^-1, 流量监测频率为1次·d^-1.TN浓度采用碱性过硫酸钾氧化-紫外分光光度法测定, 玉泉溪监测断面流量通过面积比拟法计算得到.

3.2 LOADEST模型

为了解决水质监测数据缺乏情况下的河流养分负荷量估算问题, 美国地质勘探局研发了LOADEST模型和相应的FORTRAN程序(Runkel et al., 2004).利用连续的流量和离散的水质监测数据, LOADEST模型能够自动建立逐日河流养分负荷回归统计方程, 进而可以估算日/月/季节/年等时间尺度下的养分负荷(谢润婷, 2017).该模型具有数据依赖性低, 操作简单, 结果可靠, 自动化程度高等优点.LOADEST模型有11种回归方程形式, 其中, 7参数方程表示如下：

(1)

式中, ln为自然对数, L_t为t日的养分负荷量(kg·d^-1)；Q_t为t日的流量(m³·s^-1), Q_c为中心化的流量；t为分数形式的日期, t_c为中心化的日期；α₀为常数, α₁和α₂分别指示养分负荷量与流量的关系, α₃和α₄分别指示养分负荷量与时间的关系, α₅和α₆分别指示养分负荷量与季节的关系.

LOADEST模型内置算法可以根据Akmke信息准则(AIC)和Sehwarz posterior probability准则(SPPC)自动从前9种方程形式中选择最优方程形式.模型提供了最小绝对偏差(LAD)方法、极大似然估计(MLE)和渐进极大似然估计(AMLE)3种参数率定方法, 每种方法都有其各自的适用条件(Runkel et al., 2004).

3.3 数字滤波基流负荷分割模型 3.3.1 数字滤波基流分割算法

数字滤波基流分割技术借助滤波函数将代表地表直接径流的高频信号与代表基流的低频信号分离, 实现了地表直接径流与基流的分割.数字滤波技术使用简单、结果具有可重复性和良好的客观性(Arnold et al., 1995；He et al., 2016).Nathan等(1990)首次利用该技术将基流从河川径流中分割出来, 取得了较好的结果.近年来, 随着计算机技术的发展, 源于信号分析的数字滤波技术在水文学领域得到了广泛的应用(Bosch et al., 2015；Longobardi et al., 2016；刘新有等, 2017；周星等, 2017).Eckhardt(2005)针对传统的数字滤波方程进行了改进, 引入退水常数(α)和最大基流指数(BFI_max), 提出Eckhardt′s数字滤波(ERDF), 使数字滤波基流分割技术在实际应用中更加灵活.ERDF的计算公式如下：

(2)

式中, Q_t为t时刻的径流(m³·s^-1)；q_t为t时刻的基流(m³·s^-1)；BFI_max为最大基流指数(无量纲)；α为退水常数(无量纲).

3.3.2 数字滤波基流负荷分割模型形式

河流养分负荷量取决于水体中养分浓度和河流流量.通常, 地表直接径流的流量及其携带的养分浓度受降雨量、降雨强度和人类活动等外界因素影响很大, 波动剧烈, 故地表直接径流养分出流时间序列频域下的信号较强, 因而可以看做高频信号；相较而言, 由于受到地下水库的调蓄作用, 基流流量及其携带的养分浓度相对稳定, 因此, 基流养分出流时间序列变化幅度较小, 相对平稳, 可以看作低频信号(何圣嘉等, 2017).这也意味着, 河流养分负荷量时间序列可以看作高低频信号的叠加, 换句话说, 数字滤波技术具有应用于基流负荷分割的理论依据.

本研究借鉴ERDF的形式, 数字滤波基流负荷分割模型数学形式如下所示：

(3)

式中, L_t为t时刻的河流TN负荷量(kg·d^-1)；l_t为t时刻的基流TN负荷量(kg·d^-1)；β为负荷消退常数(无量纲)；BLI_max为最大基流负荷指数(无量纲).

3.3.3 模型参数计算方法

使用数字滤波基流负荷分割模型之前, 所需解决的关键问题是计算模型参数β和BLI_max.本文建议采用如下方法计算β和BLI_max的数值.

① β通过对基流负荷消退过程的筛选和拟合得到.从经LOADEST模拟出的逐日河流TN负荷量时间序列中挑选典型基流负荷消退过程.根据如下经验公式计算峰值后直接径流持续时间, 确定某次基流负荷消退过程的起点(Linsley et al., 1949；Halford et al., 2010)：

(4)

式中, N为峰值后直接径流持续天数, A为流域集水面积.挑选出的基流负荷消退过程要满足负荷量随时间变化不断减少的特点, 并且持续天数最好大于等于5 d.在实际应用中, 由于缺乏长序列的观测数据或者有些地区降雨-产流过程频繁, 严格按照上述条件可能无法挑选出合适的基流负荷消退过程.在这种情况下, 基流负荷消退过程持续天数可以适当减少, 以便于挑选出一定数量的基流负荷消退过程, 进而计算出负荷消退常数.以L_t为因变量, 以时间为自变量, 进行指数拟合(y=be^cx), e^c即为负荷消退常数β.

② 在确定β的基础上, 借鉴基流流量分割中的反向滤波算法估算BLI_max(Collischonn et al., 2013；Miller et al., 2016).反向滤波算法原理详见文献Collischonn等(2013), 此处简要介绍BLI_max的计算过程.

地下水养分负荷的消退服从以下公式：

(5)

式中, l_t为t时刻的基流养分负荷量(kg·d^-1)；l_t-1为t-1时刻的基流养分负荷量(kg·d^-1).

用于基流养分负荷估算的反向滤波函数表达如下：

(6)

式中, l_t′为通过反向滤波得出的t时刻基流养分负荷量(kg·d^-1)；l_t-1′为通过反向滤波得出的t-1时刻基流养分负荷量(kg·d^-1).将研究时期最后一天的养分负荷量设为反向滤波算法的初始值, 可以推演出完整的基流负荷过程.BLI_max可由反向滤波基流负荷除以总负荷得出.

3.4 模型性能评价指标

采用决定系数(R²)和纳什系数(Nash-Sutcliffe model efficiency, NS)评价模型的模拟精度, 二者具体计算公式如下：

(7)

(8)

式中, n为样本数目；O和P分别代表实测数据和模拟数据的平均值；O_i和P_i分别代表第i个实测数据和模拟数据.

3.5 数据处理与图形制作

通过USGS开发的FORTRAN程序运行LOADEST模型；应用MATLAB 2014a编程实现负荷消退曲线自动挑选及反向滤波算法；其余数据处理利用Microsoft Excel 2007软件完成；采用ArcGIS 10.2和Origin Pro 9.0软件制作相关图件.

4 结果与讨论(Results and discussion) 4.1 河流总氮负荷量的LOADEST模拟

利用玉泉溪流域2010—2013年逐月TN浓度数据和相应水质监测日的流量数据, 实现LOADEST模型方程形式优选及方程参数率定.通过AIC和SPPC准则从模型常用的9种回归方程中选择最优的方程形式(Runkel et al., 2004).图 2表明, 与其他方程形式相比, 双参数方程的AIC和SPPC值最小, 表明该方程为最优的逐日TN负荷量回归方程.并且, 该回归方程残差序列相关系数(SCR)很低, 说明模拟残差之间相互独立, 方程中涉及的变量不存在序列相关性.LOADEST模型建模过程中, 参数系列a₀~a₆的取值十分重要, LOADEST模型的FORTRAN程序中内置了LAD、MLE和AMLE 3种方法用于自动计算参数系列的数值并能够分析其不确定性和统计学意义.概率曲线相关系数(PPCC)接近于1表明模拟残差符合正态分布, 因此本文使用AMLE方法进行参数率定, 各参数的率定值见表 1.从表 1可以看出, a₀和a₁的平均值分别为4.1133和0.9762, t-ratio绝对值都大于1, 表明率定后的各参数值不确定性很低；a₀和a₁与0之间存在显著差异性(p < 0.05), 具有统计学意义.由图 3可以看出, 水质监测日的TN负荷量实测值与LOADEST模型模拟值之间的R²和NS都为0.72.根据相关研究, 当NS > 0.65时, 模型的模拟结果为令人满意；当NS > 0.80时, 模型的模拟结果为相当精确(Borah et al., 2004；Shu et al., 2008).综上所述, 所建立的双参数LOADEST模型模拟效果较好, 可以有效地进行玉泉溪逐日TN负荷量估算, 形成以日为尺度的TN负荷量时间序列.

4.2 数字滤波基流负荷分割模型参数计算

根据相关公式, 确定峰值后直接径流持续时间为3 d.珊溪水源地多年平均降水天数为149 d, 降雨-产流过程频繁, 基流负荷消退过程如果根据持续天数大于或等于5 d进行筛选, 则只能从玉泉溪2010—2013年TN负荷量时间序列中选出1条负荷消退曲线.但该拟合曲线明显不符合基流负荷的消退过程, R²仅为0.62, 无法得出合理的负荷消退常数β(图 4 b).按基流负荷消退过程持续天数大于或等于4 d可以筛选出3条消退曲线, 其拟合曲线的R²分别为0.92、0.62和0.98, 计算出对应的β分别为0.74(e^-0.299)、0.23(e^-1.478)和0.85(e^-0.158)(图 4 a~4 c).因此, 本研究基流负荷消退持续天数确定为4 d.严格来讲, 实际基流消退过程中不可避免的包含一部分地表直接径流, 并且由于地表径流对地下水的补充(Eckhardt, 2008), 会造成观测到的基流负荷消退过程快于理论上的基流负荷消退过程.综合考虑上述原因及拟合曲线R²大小, 选择0.85(即基流负荷消退过程最慢)作为负荷消退常数, 用以保证所挑选的负荷消退曲线能够较好地代表流域实际的基流负荷消退过程.确定β之后, 运行反向滤波算法, 得出最大负荷指数BLI_max = 0.62(图 4 d).

4.3 基流对河流氮负荷的贡献定量

由于基流过程的复杂性、不确定性和难以监测性, 目前尚无法实现流域尺度基流流量的直接测定(徐磊磊等, 2011), 更没有直接测定基流负荷量的方法(何圣嘉等, 2017).因此本文根据分割出的基流负荷是否相对稳定及其时间序列是否符合基流非点源污染产生的特点判断计算结果的合理性和适应性.此外, 为进一步量化模型的模拟性能, 采用一种间接验证的方法对模型结果进行验证：首先通过多种基流分割方法(数字滤波法、非线性水库法和基流指数法)对玉泉溪流量进行分割并取结果的平均值；之后, 按照Schilling等(2004)提出的“纯”基流日期确定方法从水质监测日中筛选出河道水体完全由基流补给的日期, 其相应的河流负荷量近似为“实测”基流TN负荷；最后, 比较“纯”基流日期的“实测”TN负荷与模拟基流负荷的拟合程度, 利用R²和NS评估模型的性能.应用构建的数字滤波基流负荷分割模型对玉泉溪流域TN负荷量进行基流负荷分割, 结果见图 5.从图 5可以看出, 基流养分的出流时间序列相对平稳, 基流过程线较为平滑, 符合其产生的客观实际, 故玉泉溪流域基流负荷分割结果较为合理.并且, “纯”基流日的“实测”TN负荷与模拟基流负荷的拟合程度非常好, R²和NS分别达到0.83和0.80(图 6), 说明所建立的数字滤波基流负荷分割模型模拟结果相当精确, 可满足基流非点源TN负荷量定量计算的需求.

玉泉溪流域2010—2013年TN负荷量为141.21~274.68 t·a^-1, 平均208.63 t·a^-1；年基流TN负荷量为84.39~168.68 t·a^-1, 平均127.69 t·a^-1；年地表直接径流TN负荷量为56.82~106.00 t·a^-1, 平均80.95 t·a^-1(表 2).基流是玉泉溪流域氮素输出的最主要途径, 对玉泉溪年均TN负荷量的贡献率高达61.20%, 接近地表直接径流贡献率的1.5倍.本研究结果与浙江省长乐江流域(基流贡献率为61.82%)、美国Raccoon River流域(基流贡献率约为66%)、韩国Daejeon市地表水(基流贡献率约为60%)和美国Chesapeake湾流域(基流贡献率为58%~73%)等地区研究结果相似(Schilling et al., 2004；Kim et al., 2010；Miller et al., 2016；何圣嘉等, 2017).这一现象主要由以下原因引起：①基流是玉泉溪主要的水源, 基流对总径流贡献率达50%以上, 尤其在枯水季节, 河道水体通常都由基流补给；②流域内林地面积约占流域总面积的65%, 流域范围内广泛分布的林地能够拦截、降解和吸收地表直接径流携带的养分, 可以有效降低地表直接径流TN负荷量.由基流负荷分割结果可见, 当前玉泉溪氮负荷主要来源于该流域基流养分的输出, 故而在库区非点源污染的治理工作中, 只有充分认识到地下水对地表水污染的重要性, 才能制定更加科学合理的流域养分管理方案.

4.4 玉泉溪流域氮污染控制对策

明确基流输出对河流氮污染的贡献能够指导相关部门从短期和长期的角度制定有针对性的养分控制措施, 合理分配有限的资源用以提高水污染治理的效率.许多地区非点源污染理论研究和控制实践表明：在基流输出成为当前地表水氮污染的重要甚至主要来源的情况下, 通过削减氮素源头输入控制水体氮污染的成效有可能要经历较长的时间才能呈现.例如, 有研究显示, 丹麦Kastbjerg和Villestrup流域从20世纪80年代初期实施了严格的氮素削减措施, 农业氮素排放已经显著减少, 但是此后15年河流水体氮素浓度却仍然呈现上升趋势(Vervloet et al., 2018).Van Meter等(2018)研究结果表明由于地下水系统遗留硝态氮的输出, 美国Mississippi河流域即使从2017年开始通过一系列措施将农业活动中的氮素利用率提高到100%(流域氮素盈余为0), 也需要经过约30年的时间才能达到削减60%硝态氮输出负荷的预期目标.针对水源地非点源污染问题, 从长期的角度而言, 通过退耕还林、测土配方施肥、提高水肥管理效率和保护性耕作等“源头控制”措施, 从源头上减少进入环境的养分是非点源污染防控的根本.从短期来看, 由于养分通过地表直接径流和基流入河的时间和途径不同, 需要对二者加以区分, 以便有针对性地采取相应的措施迅速改善河道水质状况.例如, 地表直接径流污染贡献率较高的流域, 意味着水体氮素来自滞后性污染源的比例低, 因而河流水质对源控制策略的响应较快, 可实施“养分源头控制”.并且, 生态沟渠、植被过滤带、生态护岸等“陆地过程拦截”措施能够减少水土流失、对地表直接径流具有较好的净化效果, 进而削减地表直接径流带入受纳水体的养分.在此基础上, 辅以“河道过程拦截”措施能够进一步提高水污染治理的成效.对于基流氮素输出较高的流域, 采用“养分源头控制”和“陆地过程拦截”措施很难在短期内实现水体氮污染的改善.在目前还没有特别有效的办法用于地下水氮污染治理的情况下, 通过“河道过程拦截”措施强化河道对氮素的净化功能, 从而迅速减低氮素浓度, 应该作为非点源氮污染前期治理的优先选择措施.

玉泉溪流域基流输出贡献了河流TN负荷量的60%以上, 这意味着该流域河流氮污染有大部分来自于地下水中累积的养分.因此, 在该流域非点源污染控制实践中, 从短期来讲, 河道过程拦截措施相对于源头控制更能迅速地减少河流养分入库负荷量, 从而降低珊溪水库水体富营养化的风险, 保障饮用水供应安全.对于不同类型的河道, 应该因地制宜地构建河道过程拦截措施, 充分发挥生物对水体中氮素的吸收降解作用, 提高水体的自净能力：考虑到玉泉溪入库河口区水面宽广, 水流速度慢, 因此可以采用人工浮床和前置库等措施, 在玉泉溪入库之前拦截削减水体中氮素, 减少入库氮通量；下游河道(尤其是入库前河道)水面较广, 水位较浅, 可以采用生物修复(沉水植物、挺水植物)、人工湿地和人工浮床等措施强化河道水质净化功能；中上游有些河道狭窄, 可以通过加设橡胶坝和拓宽河道等工程措施增加水体停留时间, 进而一定程度上提高河道的净化功能.本研究结果也从一定程度上解释了尽管库区环境综合整治工作不断推进, 但水体氮/磷污染没有得到根本性好转的原因可能在于流域地下水中累积的养分库通过基流过程的不断缓慢输出(Sanford et al., 2013；Vervloet et al., 2018).需要注意的是, 由于地下水更新过程持续时间较长, 该流域水体质量的根本性改善将会是一个长期的过程, 需要正确引导政府和民众对于水污染治理工作的预期, 以保持对水源地水质改善的信心和投入.

5 结论(Conclusions)

1) 本文充分考虑了水源地基础数据缺乏的现状, 基于LOADEST模型和数字滤波算法创建了地表直接径流与基流养分负荷分割的方法体系.经过实例验证, 所提出的水源地基流氮素定量方法模拟得到的基流养分出流时间序列符合其产生的客观实际, 并且评价指标R²和NS分别达到0.83和0.80, 为水源地基流养分负荷的估算提供了可靠工具.所建立的方法具有以下明显的优点：可操作性强, 建模简单方便, 易于推广使用；大大的减少了对数据的需求, 只需要大多数地区都能很容易获取的常规水质和流量数据就可以实现基流负荷量的计算；数字滤波基流负荷分割模型中滤波参数是根据研究区养分负荷时间序列自身特点进行计算, 避免了参数选择的主观性, 结果更加符合不同研究区的客观情况.

2) 玉泉溪流域基流负荷分割结果表明基流氮素输出是造成玉泉溪TN污染的主要原因, 2010—2013年均基流对河流TN负荷量的贡献率为59.76%~61.97%, 平均61.20%.因此, 在玉泉溪流域饮用水源保护工作中, 应该充分考虑流域养分以基流形式的输出, 制定更加科学高效的非点源污染控制方案.