洪泽湖入湖河流对湖区水质的响应关系

刘 超， 胡智华， 姚天启， 李 港

(1.重庆交通大学河海学院， 重庆 400074； 2.南京水利科学研究院生态环境研究所， 江苏 南京 210029；3.贵州民族大学生态环境工程学院， 贵州 贵阳 550025)

洪泽湖流域地处经济高速发展的苏北平原，流域内工农业以及水产养殖业发展迅速，为区域经济发展做出了强有力的贡献，由此衍生出较为严重的环境问题，如地表水环境受到强烈的人为活动干扰，氮磷质量浓度远高于目标水质标准。洪泽湖属于典型过水型湖泊，入湖河流是湖区水资源的主要补给方式，因此改善湖区水质最关键手段是控制外源输入。从1990—2018年的时间过程中，洪泽湖水质整体呈现出V类到劣V类再到V类的变化趋势，主要水质超标指标为TN和TP[1]。研究表明2000年以来洪泽湖实施的各种治理和管控措施对水质改善具有一定的成效，但湖区整体水质仍然未达到III类目标水质标准，且空间差异仍然十分明显。因此，进一步开展降低洪泽湖水质空间差异和改善洪泽湖整体水质的研究具有重要意义。

洪泽湖流域涉及淮安市和宿迁市的6个区县，各个区县经济发展水平以及产业侧重点各不相同，导致污染物的来源存在较大差异。同时受不同入湖河流水量以及湖区自身地形复杂致使湖泊各个区域水动力条件有着较大差异，各入湖河流对洪泽湖不同区域水质的影响程度不一。目前关于洪泽湖水质改善的研究相对比较宏观，如徐嘉兴等[2]研究了洪泽湖流域土地利用组成，发现加大集中生活废水处理力度以及减少水产养殖面积，可降低入湖河流污染负荷。阎敏等[3]研究表明农业面源是入湖河流污染负荷主要来源，控制耕地面积有利于湖区水质的改善。已有研究都相对缺乏以洪泽湖水质作为目标进行精细化的探索，只有识别入湖源与湖区水质在时空上的定量响应关系，才能更加准确实施控源措施以及提升水质的改善效率。因此，准确识别各入湖河流与湖区水质在时间和空间上的定量响应关系亟待解决。

本研究通过环境流体动力学模型(Environmental fluid dynamics code, EFDC)建立洪泽湖水动力水质模型，将概化的7条主要入湖河流作为边界，选取洪泽湖主要的水质超标指标TN和TP作为模拟输出目标，探究不同入湖源对湖区水质的定量响应关系，优化洪泽湖外源输入控制措施，为提高洪泽湖水质改善效率提供科学的理论依据[4]。

1 材料与方法

1.1 研究区域

洪泽湖是我国第四大淡水湖，属于典型过水型湖泊，是淮河流域最大的防洪供水调蓄湖泊，也是南水北调东线工程的重要水源地，在调节气候和保护生物多样性等方向发挥着极其重要的作用。洪泽湖气候四季分明，多年平均降水量为925.5 mm，每年6—9月降水量占全年总降水量的65.5%，冬季降水量仅占全年7%～8%。洪泽湖流域保护面积达到3 704 km2，湖泊保护范围线长1 006 km。主要入湖河流有淮河、成子河、濉河、老濉河、新汴河、怀洪新河、徐洪河，淮河是洪泽湖主要的供水河道，出湖河流主要有淮河入江河道和苏北灌溉总渠等。

1.2 水质数据

江苏省水文水资源勘测局在洪泽湖水质监测站点共15个，#1、#2、#3、#5监测站点位于成子湖区，#4、#6、#7、#8监测站点位于湖心区，#12、#13、#14、#15监测站点位于溧河洼区，#9、#10、#11监测站点位于入江过水通道区，各个监测站点的监测频次为12次/a(每月1次)。本研究数据时间序列为2016—2018年，监测指标为DO、水位、TN和TP，DO采用水质监测仪直接测得，水位为当地水位站实时监测数据。TN和 TP为现场分层(表、中、底)取混合水样10L冷藏运回实验室检测，TN采用碱性过硫酸钾消解紫外分光光度法(CB 11894—1989)，TP采用钼酸铵分光光度法(GB 11893—1989)。

1.3 情景模拟方案设计

目前关于洪泽湖水质变化特征的研究相对较多，已有研究表明洪泽湖水质在时间和空间上都有着显著的差异。利用2010—2019年洪泽湖长序列水质监测数据分析得到，TN和TP年平均质量浓度变化范围为1.39～1.86mg/L和0.08～0.171mg/L，变化幅度分别为33.81%和26.25%。为准确探究入湖源在时间和空间上对湖区水质的定量响应关系，避免水质在时间上的波动以及水质模型自身存在的误差导致模拟结果变化不显著的影响，对入湖源TN和TP质量浓度以50%削减率作为标准情景工况，采用一次一个变量的单一控制变量法对7个情景方案进行情景模拟，分析计算不同入湖源对湖区不同监测点位TN和TP质量浓度的改善能力。

2 水动力水质模型

2.1 模型原理

EFDC模型水质模块[5]中涉及氮的变量共有5个，分别为难溶性颗粒态有机氮(RPON)、活性颗粒态有机氮(LPON)、溶解态有机氮(DON)、氨氮(NH4)和硝态氮(包括硝酸盐氮：NO2；亚硝酸盐氮：NO3)。涉及磷的变量有4个，分别为难溶性颗粒态有机磷(RPOP)、活性颗粒态有机磷(LPOP)、溶解态有机磷(DOP)和TP酸盐(PO4)。模型中各水质变量彼此之间具有如下关系结构。

本研究运用的水动力水质模型是在Sigma坐标下基于垂向静水压力的假设，并采用布西内斯克方程(Boussinesq equation)近似得到如下的动量方程和连续性方程。

x方向的动量方程为：

(1)

y方向的动量方程：

(2)

z方向的动量方程:

(3)

连续性方程：

(4)

(5)

密度方程：

ρ=ρ(p,S,T,C)

(6)

(7)

P=myHuQ=mxHv

(8)

(9)

式中：(x，y)为水平方向的曲线-正交坐标；z为垂向σ坐标；(u，v)为(x，y)方向的水平速度分量，m/s；H为总水深，m；mx，my分别为x，y坐标变化系数，在笛卡尔坐标下，变化系数等于1；Patm为大气压强，Pa；p为参考密度ρ0下的附加静水压；b为浮力；f为科里奥利力系数，涵盖网格曲率加速度；AH为水平动量扩散系数，m2/s；Av为垂向紊动粘性系数，m2/s；cp为植被阻力系数；Dp为与每单位水平面积的流量相交的投影植被区域；Su和Sv为(x，y)方向动量方程的源/汇项，m2/s2；Sh为质量守恒方程的源/汇项，m3/s；S为盐度，ng/L；T为温度，℃；C为总悬浮无机颗粒质量浓度，g/m3；U和V为(x，y)方向的质量浓度平均速度分量，m/s；P和Q为(x，y)方向的质量通量分量，m2/s，w为σ坐标下的垂向速率，m/s；w*为z坐标下的垂向速率，m/s。

2.2 模型构建

本研究所建立的EFDC水动力水质模型是耦合了水动力和水质两个模块。洪泽湖岸线提取，利用正交网格将整个洪泽湖划分成9 749个矩形网格单元，每个网格单元边长为400 m×400 m，垂向采用Sigma坐标系，平均分为3层，各层占水深比例为0.333，为了适应水位波动带来的影响，设置临界干水深为0.07 m。边界条件有气象边界条件、水动力边界条件以及水质边界条件。气象边界条件主要包括降水、气压、气温、太阳辐射、湿度、蒸发量以及云覆盖系数，气象资料均来自洪泽湖周围气象站，获取方式为中国气象数据网。洪泽湖周边河道复杂多样，省略一些较小且数据不完善的河流，最后将主要入湖河流概化为7条，流量边界采用2018年1月1日到2018年12月31日的逐日流量数据，数据来源为水文统计年鉴。选用洪泽湖水质主要超标项目TN和TP作为模型输出目标，水质数据包含入湖河流和湖区15个监测站点2016—2018年月尺度数据。初始水位和水温均用2018年1月1日实测值。模型初始计算日期为2018年1月1日，运行周期为364 d，时间步长为20 s。

2.3 模型参数敏感性分析

根据输出目标所涉及的模块，通过敏感性分析选择敏感性较高的参数，对所选择的参数进行敏感性分析，确定其最佳的取值。通过已有文献确定模型参数的最小值与最大值，通过LHS(Latin hypercube sampling)方法抽样分布函数，基于参数分布范围利用LHS生成20组随机参数组合，针对LHS生成的每一组参数，运用模型进行模拟，输出目标结果，评估参数不确定性对输出目标结果的不确定影响，最后采用标准秩回归法进行参数敏感性分析，确定各个参数的最优取值，如表1所示[6]。

表1 模型校正后主要参数取值

2.4 模型率定和验证

为确保水动力水质模型的可靠度，需进一步对模型进行率定和验证，模型验证的目标变量为水位和DO，模型率定数据采用2016年1月1日至2016年12月31日实测数据，模型验证数据采用2017年1月1日至2017年12月31日实测数据。

根据湖区3个水位站(高良涧站、蒋坝站、老子山站)实测水位数据与模拟水位数据对比分析可以发现模拟值与实测值较为接近，变化趋势比较同步，模拟与实测数据决定系数(R2)均高于0.918，这说明模型在水位模拟上能够较精确的反映真实水位情况，这为TN和TP等水质指标的模拟提供和可靠地支撑。利用湖区3个水质监测站点(#5、#14、#8)的DO实测数据与模拟值对比分析发现，模拟值与实测值的变化趋势以及差值大小基本在可接受范围之内，说明建立的湖区水质模型能够较好的反映出湖区水质的变化和分布情况，其模拟值与实测值决定系数(R2)最小为0.8246，最大值为0.8722。水动力R2略高于水质的R2，这表明水质在模拟过程中的输出结果相对于实测值的误差比水动力模拟输出结果相对于实测值的误差要稍大一些，出现这种情况主要是因为水体中水环境因子变化复杂、彼此之间相互影响，加上水质的模拟结果是建立在复杂的水动力基础之上，受到水动力的影响。

3 结果与讨论

3.1 入湖源与湖泊整体水质的响应关系

洪泽湖各入湖源总入湖水量、总污染负荷以及在时间进入湖区的水量分配皆有着较大差异，6—8月为丰水期，进入湖区水量最大，10月为平水期，12月为枯水期，进入湖区水量相对较小[7]。入湖水量最大为淮河入湖河流(301.4 × 108m3)，TN年平均质量质量浓度最高为新汴河(2.82 mg/L)入湖河流，TP年平均质量质量浓度最高为老濉河(0.183 mg/L)入湖河流。

图1 水位和DO率定结果

为探究7个入湖源对湖区水体的影响关系，以代表性、全局性为原则在湖区布设了15个监测站点。本研究以单一变量的研究方法进行了8种情景工况(1组未削减和7组削减情景工况)的模拟，得到不同情景工况下湖区水体TN和TP的质量质量浓度的变化过程，并以15个监测站点平均质量质量浓度来表示湖区水体营养盐的整体质量质量浓度水平，从而分析入湖源对湖区整体水质的影响能力。

削减模拟结果与未削减的基础模拟结果如图2所示。从时间上看，模拟前20 d各削减情景工况与未削减基础情景工况的模拟结果未出现显著差异，这主要是因为模型运行初期模型尚未稳定，故而受到模型初始条件的影响。初始条件是根据模拟年份第一天的实测数据进行设定，只有在模型运行稳定，初始条件的影响消失过后计算出的结果才可靠。不同入湖源进行相同比例削减模拟结果显示，淮河对湖区TN和TP年平均质量浓度改善效果皆为最好，主要是因为淮河入湖流量大，流速快、水动力条件好，因此污染物迁移能力强。

3.2 源边界影响特征的时空间差异性

对于相同的入湖源，在不同的时间节点上，入湖源对水质的改善效果也有差异，反之，不同的入湖源在相同的时间节点上亦如此。为探究不同入湖源在各自对应的不同时间上对TN和TP质量浓度的改善差异，分别选取了位于成子湖、溧河洼、湖心区、入江过水通道的4个监测站点对应第210天(丰水期)、第300天(平水期)、第350天(枯水期)进行分析。结果如图3～ 5所示。

在相同时间节点上各入湖源对不同监测站点TN何TP的改善效果不同。如图3中R3入湖源对#2号监测站点TN改善率占比为40%，对#14号监测站点TN改善率占比为4%，对#6号TN改善率占比为14%，对#9号TN改善率占比为21%。R3入湖源对#2、#14、#6、#9号TP改善率占比分别为65%、3%、5%、11%。同时可以发现入湖源对监测站点水质改善率主要表现为近高远低的响应特征。但也有不同的响应规律，如图3中R1入湖源对#6和#9号监测站点的改善率占比分别为65%和3%，表现出近低远高的响应特征，这主要是因为丰水期淮河入湖流量剧增、污染负荷增加，流速增大，在剧烈的水动力条件下迁移方向较为径直，导致出现对较远径直站点影响更大[8]。这意味着监测站点水质改善效果主要受到入湖源水量、流速、污染负荷量以及监测站点与入湖源之间距离的综合影响。

在不同时间节点上，入湖源对同一个监测站点改善效果不同。如图3～5中R1入湖源对#14号监测站点TN的结果显示，在第210天时改善率占比为68%，第300天时改善率占比下降到55%，第350天时改善率占比进一步下降到42%。图3～5中R1入湖源对#14号监测站点TP的结果显示，在第210天时改善率占比为59%，第300天改善率占比为35%，第350天改善率占比为58%。这主要是因为在不同水情条件下，入湖源的水量和携带的负荷总量的动态变化过程明显[9-10]。这表明对同一个监测站点而言，充分考虑入湖源的削减时间对于提高水质改善率有着正向效应。对同一监测站点在不同时间入湖源对其改善率不同。如图3～5中就#9号监测站点TN而言，在第210天时改善率最大为R5入湖源(26%)，在第300天时改善率最大为15%，在第350天时改善率最大为R1(44%)。这意味着，提高同一监测站点的水质改善率，在关注时间的同时也要对入湖源进行有针对性的筛选。

图3 TN和TP模拟210 d 4个监测站点解释结果

图4 TN和TP模拟3 000 d 4个监测站点解释结果

图5 TN和TP模拟350 d 4个监测站点解释结果

4 结 论

淮河入湖源对洪泽湖水质指标TN和TP平均质量浓度的改善效果最好，改善率分别为16.2%，22.95%；徐洪河入湖源对最大值质量浓度改善效果最好，改善率分别为17.34%，24.08%。当以平均水质作为目标时，对淮河进行控源效果更明显，当以降低水质峰值质量浓度为目标时，因对徐洪河和成子河入湖源实施控源措施。时间差异性，不同时间节点入湖源对同一个监测站点改善效果不同。对于成子湖在丰平水期对徐洪河进行控源对水质就能有较好的改善效果，在枯水期需要对淮河与徐洪河联合控源才能有较好的改善效果。空间差异性表明，入湖源对监测站点水质改善率主要表现为近高远低的响应特征。环湖边界水质受就近入湖源影响较大，向湖心区移动受淮河入湖源影响逐渐增强。入湖河流对洪泽湖水环境健康有着极其重要的影响，在以改善湖区水质为目标的管理中，应当综合考虑各个入湖河流总水量、流速、水量的时间分配、污染物质量浓度等特征，根据入湖河流与湖区水质在时空上的定量响应关系做出精细化的控源措施。