洪泽湖换水能力的时空分布特征

薛联青，沈海岑，张 敏，汪 露

(1.河海大学水文水资源学院，江苏 南京 210098; 2.皖江工学院水利工程学院，安徽 马鞍山 243000；3.江苏省洪泽湖管理委员会办公室，江苏 淮安 223100)

洪泽湖是淮河流域重要的调蓄湖泊，淮河下游地区重要的水源地，同时也是南水北调东线工程的主要调节湖泊[1-2]。随着区域人类活动加剧，洪泽湖面临着自由水面面积缩减、局部富营养化、植被覆盖度降低等多重环境问题[3-8]。湖泊的换水周期决定着输移容量和水环境容量，是判断湖泊污染物输移能力的重要指标[9]。风场和吞吐流对浅水湖泊的水体交换过程有重要影响[10]，洪泽湖的湖面宽阔，风场和吞吐流共同塑造其流场形态[11]，从而对洪泽湖污染物的扩散、迁移产生影响。因此探究不同风场和不同季节来水条件下洪泽湖换水周期的时空分布，有利于识别水动力因素对污染物运移的影响[12]。

国内外众多学者对湖库、海湾等水体换水能力开展了研究。Cucco等[13]构建了威尼斯潟湖二维水动力模型，通过潟湖内释放的被动示踪剂的残差函数来确定水体停留时间；Bocaniov等[14]利用三维水动力-富营养化模型计算了3种水体传输时间，研究发现营养盐损失率与所有运输时间尺度指标均呈统计学显著相关；李云良等[15]基于对流扩散理论计算不同季节鄱阳湖的换水周期和传输时间，研究发现季节水情对换水周期有较大影响；胥瑞晨等[16]研究得出水龄和水体停留时间存在63%的比例关系，并计算出太湖在不同风场条件下不同方位湖区的水体交换率和半交换周期；陶磊等[17]基于欧拉观点建立了渤海湾水龄计算模型，分析了潮汐和季风两个因素对渤海湾水体交换的影响。目前，对洪泽湖水体交换能力的研究仅停留在定量计算阶段[18]，缺少对水体交换驱动机理的研究。本文基于MIKE21建立洪泽湖二维水动力模型，模拟计算洪泽湖不同湖区的换水周期，考虑洪泽湖的流场呈现风生流和吞吐流双重特点，分析不同季节和不同风场条件下换水周期的空间变化特征，从流场的角度分析水体交换的驱动机理，可为洪泽湖的水环境治理提供科学依据。

1 研究区概况

洪泽湖的多年平均水位(蒋坝站)为12.55 m，平均水深1.9 m，最大水深4.5 m。多年平均入湖水量为342亿m3，入湖河流主要分布在湖泊西岸和南岸，包括淮河干流、怀洪新河、徐洪河、濉河、新汴河、老濉河和老汴河(图1)，其中，淮河干流每年入湖水量超过70%。多年平均出湖水量为313亿m3，出湖河道为二河、入江水道和苏北灌溉总渠，其中，入江水道为主要出湖河道，控制工程为三河闸，出湖水量占总量的60%～70%。洪泽湖的不同湖区呈现不同的流场形态，流场形态对大型过水性湖泊的污染物运移方式有重要影响[19]，不同湖区的水质状况差别较大，结合湖泊特点将洪泽湖划分为成子湖区、二河区、临淮区、老子山区、三河区、中心湖区和溧河洼区7个湖区，如图1所示。

图1 洪泽湖分区及水系分布

2 研究方法

2.1 水动力模型原理

MIKE21水动力模型基于Navier-Stokes方程，并遵从Boussinesq假定。二维非恒定流浅水方程为

(1)

(2)

(3)

其中

2.2 水动力模型建立及验证

采用洪泽湖2019年90 m×90 m的DEM数据，结合实测水深资料，在MIKE模型中进行修正，得到湖底高程数据，利用Mesh Generator工具将湖区划分为19 475个非结构三角形网格，共计10 888个节点。将洪泽湖的主要入湖口和出湖口作为开边界，综合考虑洪泽湖水系特征，入湖河道选择淮河干流(包括池河)、怀洪新河、新汴河、濉河和汴河，出湖河道选择三河、二河、苏北灌溉总渠和徐洪河。上边界(入湖边界)输入实测流量数据，下边界(出湖边界)输入实测水位数据。采用2016年的实测水位数据对模型进行验证，模拟时段为2016年1月1日8：00至2016年12月31日8：00。初始水位设为2015年12月31日蒋坝站水位13.12 m，初始流速设为0 m3/s，最小时间步长设为1 s。根据洪泽湖实际情况，设定干水深为0.006 m，湿水深为0.1 m，淹没水深为0.05 m。主要参数设置为：曼宁系数45 m1/3/s、CFL数0.8、涡黏系数0.28。验证结果见图2，计算求得尚咀、老子山、高良涧闸和蒋坝4个验证点的Nash系数分别为0.988 6、0.990 2、0.995 3和0.988 4，满足模拟的精度要求。

图2 模型验证结果Fig.2 Results of model validation

2.3 换水周期计算原理

很多学者从不同角度定义湖库的水体交换能力，比如水龄、示踪剂传输时间、换水周期等[20-22]。湖泊换水周期是评价水体中由对流扩散导致的物质交换速率的常用概念，本文根据湖泊特点，选择湖泊换水周期来定量研究洪泽湖的水体交换能力。在湖泊中投入不可降解的示踪剂，假设示踪剂浓度变化呈现指数规律衰减[23]，则湖泊平均换水周期计算公式为

(4)

其中

式中：τ为湖体平均换水周期；rt为t时刻示踪剂质量浓度与初始时刻质量浓度之比；ρt为t时刻示踪剂质量浓度；ρ′0为初始时刻示踪剂质量浓度。

由于现实中物质浓度不可能降为0，当t=τ时，示踪剂浓度降为初始浓度的e-1或37%[24]，定义示踪剂浓度衰减至初始浓度的37%所需的时间为换水周期。计算换水周期时，将初始示踪剂空间质量浓度场设为1 mg/L，入湖口质量浓度设为0 mg/L，出湖口设为自由出流，通过编写Python3.9代码读取每个非结构网格点质量浓度降解为0.37 mg/L的时刻，用空间连接功能导入到ArcGIS10.5中，生成湖泊换水周期的空间分布，并通过分区进行区域空间统计，计算出各个湖区的平均换水周期。

2.4 模拟情景设置

水文和风场条件是研究浅水湖泊水体交换过程的基础[25]。为研究不同季节、风向对洪泽湖水体交换能力以及污染物运移的影响，设置几种数值模拟情景(表1)，模拟时间均为365 d，即1个水文年。情景1模拟2013—2019年平均来水、气象条件下的换水周期，将示踪剂(初始浓度场)的投放时间设定为1月1日、4月1日、7月1日和10月1日，表征洪泽湖春、夏、秋和冬季的换水能力[17]。情景2设定实际风场、无风和典型主导风场为1.88 m/s的东南风与东北风，模拟不同风场条件对湖泊各区域换水周期的影响。

表1 模拟情景设置

3 结果与分析

3.1 换水周期时空分布规律

图3为不同季节换水周期空间分布，图4为各湖区不同季节的换水周期。由图3、图4可见，洪泽湖换水周期空间异质性明显，整体呈现出南部湖区水体更新时间快、北部湖区水体更新时间慢的特点，换水周期从南到北的梯度变化，主要与洪泽湖的水文情势密切相关，淮河干流入流超过入湖流量的70%，且从南部老子山处入湖，故淮河上游的入流增强了南部湖区的水体交换能力。

图3 不同季节换水周期的空间分布

图4 各湖区不同季节的换水周期

冬、春两季全湖平均换水周期较长，分别为75 d和60 d，秋、夏两季的全湖平均换水周期较短，分别为49 d和31 d。冬、春两季的换水周期空间分布相似，淮河入湖口(老子山区)湖泊换水能力较强，换水周期约为10 d，呈现从淮河入湖区到二河、三河出湖区递增的特点，三河区和中心湖区的换水周期均值为45 d左右，二河区的均值为83 d左右，出湖口局部水流交换较快，在10 d以内。成子湖区的换水能力最弱，湖水流动缓慢，均值为157 d左右，成子湖区最北部的换水周期最长，达250 d。西部湖区溧河洼区和临淮区在冬季的换水周期均值分别为80 d和46 d，春季较冬季换水能力增强，换水周期分别缩短了31.25%和37.5%。

夏季属于丰水期，淮河干流入湖流量占全年的40%左右，上游大量来水导致全湖水体更新能力显著增强，换水周期均值为31 d，比冬季和春季减少了近50%，南部湖区均值小于10 d。其中溧河洼区、中心湖区、临淮区和三河区水体更新能力较春季改善显著，换水周期均从30～60 d减少为10 d左右，东部沿岸为10～20 d；二河区换水能力较弱，为30 d左右，但较春季和冬季减少了63%；成子湖区的换水周期均值为103 d，湖区最北端的换水周期长达239 d，但总体来说，均值较其他季节减少了34.8%，其中，徐洪河入湖口的换水周期小于10 d，成子湖区南部沿岸小于90 d，中心区域为90～150 d。

秋季为退水期，全湖的换水周期均值为49 d，较夏季增加了18 d。老子山和溧河洼湖区的换水周期总体仍在10 d以内，三河区、临淮区和中心湖区总体为10～20 d，二河区为46 d，较夏季增加了16 d。成子湖区换水周期分布的空间异质性明显，南部徐洪河入湖口换水周期在10 d以内，该区域面积较夏季增大，最北部的换水周期较其他季节明显增加，最长达298 d，超过200 d的湖泊面积较其他季节显著增加，这可能是由于南水北调东线工程运行之后，徐洪河作为南水北调东线的复线河道，承担抽调洪泽湖水向北送的重要任务，上半年在抽调水的过程中，徐洪河口都是水流出湖状态，根据多年平均流量显示，7月开始徐洪河入湖口水流由出湖转变为入湖，且10月后入湖流量增加，达100 m3/s，由于徐洪河湖区整体流速较缓，且风生流在成子湖区有显著作用，常出现风生环流，因此夏季、秋季在徐洪河入湖口形成换水周期小于10 d的条带，成子湖北部区域受环流影响较大，换水周期仍较长，超过200 d。

3.2 换水周期对风场的响应

风场对换水周期的空间分布会产生一定的影响[26-27]，图5为不同风场作用下的换水周期，图6为各湖区不同风场作用下的换水周期。由图5、图6可见，换水周期总体呈现从东北到西南递减的特征，对比实际风场和无风时洪泽湖换水周期的空间变化，可验证风场是影响洪泽湖污染物运移的重要驱动力。对比1.88 m2/s的东南风和东北风作用下的换水周期分布，可以看出不同的风向使得换水周期发生相应变化。

图5 不同风场作用下的换水周期

图6 各湖区不同风场作用下的换水周期Fig.6 Water exchange cycle under different wind fields in each lake area

总体来看，实际风场作用下，湖泊换水周期的均值为72 d，比无风时减少了16 d，说明实际风场对洪泽湖的换水能力起到正向作用，风生流可以促进污染物的运移扩散。风场对洪泽湖区换水周期空间分布影响最大的是溧河洼区，在实际风向作用下，其换水周期均值为92 d，比无风时减少了89 d，说明实际风场加速了溧河洼区的水体流动，增强了该区域污染物的稀释能力。临淮区的平均换水周期为55 d，减少了21 d，临淮区和溧河洼区均位于湖西，可以看出在无风条件即仅受吞吐流作用时，湖泊水体整体向东南、东北方向流动。中心湖区、二河区和成子湖区的平均换水周期在实际风场的作用下分别为36 d、67 d和204 d，仅分别减少了5 d、5 d和7 d，而三河区和老子山区为40 d和10 d，增加了14 d和2 d，说明东北部湖区的换水能力受风场的改善效果较小，主要进湖、出湖区域的换水能力略微减弱，其中，老子山区受风场的影响最小，说明该区域的水体更新主要受吞吐流的影响。

在东南风的作用下，湖泊换水周期的均值为73 d，比无风时减少了15 d，说明东南风对洪泽湖的换水能力起到正向作用，东南风可以促进污染物的运移扩散，从而加强了污染物的稀释能力。受东南风影响最大的区域是溧河洼区，其换水周期均值为58 d，比无风时减少了123 d，说明东南风加速了溧河洼区的水体流动，增强了该区域污染物的稀释能力，同属西部湖区的临淮区为48 d，减少了28 d，换水能力也明显增强。东南风也显著改善了成子湖区的平均换水能力，其换水周期为158 d，较无风时减少了53 d，最高值为193 d，较无风状态减少了47.1%。成子湖区流速缓慢，污染物容易在此聚集，受吞吐流影响较弱，所以湖流主要受风场影响，东南风加速了成子湖区的水体更新。东南风对其他湖区的换水能力产生了减弱的作用，其中，受影响最大的是二河区，其换水周期均值为121 d，较无风时增加了49 d。其次是三河区，换水周期均值从26 d增加到52 d，换水能力减弱了50%，中心湖区的换水周期也增加了17 d，由41 d增加到58 d。同时，老子山区的均值也由8 d增加到16 d，老子山区的淮河入湖区域的换水周期仍然为1～10 d，而其东部靠近三河区的区域换水周期最高值为58 d，导致均值增加了8 d。东南部湖区换水能力普遍减弱，可能是由于东南风的影响，淮河干流由老子山进入湖区后向三河区流动的水流受到阻碍，水体向东北方向流动的趋势增强，加强了东北部水体污染物浓度的稀释能力。

在东北风的作用下，湖泊换水周期的均值为106 d，比无风时增加了18 d，说明东北风对洪泽湖的换水能力起到抑制作用，东北风减缓了污染物的运移扩散速率，从而减弱了污染物的稀释能力，同时，各个湖区换水周期的整体梯度性分布特征明显，呈现东北高值，西南低值的局面。受东北风影响最大的是临淮区，其换水周期均值为145 d，比无风时增加了69 d，湖泊在无风状态下，受吞吐流的影响，淮河干流入湖流向临淮区，东北风起到抑制作用，减弱了临淮区对污染物的稀释能力。二河区的换水能力也显著减弱，其平均换水周期由76 d增加到145 d；中心湖区从41 d增加到69 d；三河区从26 d增加到43 d。东北风仅对成子湖区的水体更新起到正向作用，但对成子湖区平均换水能力的改善效果较小，成子湖区的换水周期均值为205 d，较无风时仅减少了2.8%；对换水周期最高值所在点位的水体改善效果较好，无风时，成子湖区最北边的水体交换周期接近365 d，东北风作用下最高值为246 d，减少了32.6%。东北风也降低了溧河洼区的换水周期，溧河洼区的换水周期为152 d，较无风时减少了16%。东北风对老子山区的影响较小，其平均换水周期由8 d增加到11 d，进一步验证了老子山湖区的换水能力主要受吞吐流影响。

3.3 典型湖区流场特征

风场、吞吐流以及湖泊形态影响洪泽湖流场的结构，流场结构对大型浅水湖泊的污染物运移过程和水体更新能力都会产生重要影响[28-29]。选取换水周期最短的老子山区以及换水周期最长的成子湖区为典型湖区，从流场的角度分析洪泽湖换水周期时空变化的机理。

老子山区的流场主要受淮河干流的入流影响，由图7可见，较枯月份淮河干流流量约为400 m3/s，老子山区的流速为0.024～0.048 m/s，较丰月份淮河干流流量约为1 500 m3/s，老子山区的流速为0.09～0.195 m/s。在枯水期，来水较少，流速较慢，老子山区水质变化受外界输入的影响减弱，受水体自净的影响增强，而在丰水期，流速提高了2倍，外界输入对老子山区水质变化的贡献度提高至90%。在来水较丰时期，淮河干流入湖流量大，导致老子山区水体流速加快，从而加快污染物的稀释，提高湖泊的换水能力。淮河干流的水质状况是影响老子山区水质的重要因素，但由于水体流速快，水质变化受底泥释放的影响较小，营养盐不易聚集，局部富营养化的风险较低。

(a) 2018年1月17日 (b) 2018年7月30日

由前文计算结果可知，成子湖区的换水周期大部分区域均超过150 d，模拟得到丰、枯来水条件下成子湖区的流场如图8所示。由图8可见，来水较少时，成子湖区大部分区域的流速为0.005～0.03 m/s，流速缓慢，受吞吐流影响较弱，主要受风场影响，局部呈现逆向环流；来水较丰时期，成子湖区流速略微加快，大部分区域流速为0.012～0.04 m/s，仍然存在局部逆向环流。说明无论是丰水期还是枯水期，成子湖区受吞吐流的影响较小，流速缓慢，造成成子湖区的换水周期超过100 d，最北部湖区超过200 d。由于水流交换过缓，可以推测出外界输入污染负荷对成子湖区水质变化的影响较小；同时，流速缓慢的水动力条件下，水体中的氮磷营养盐主要呈现溶解态，浓度高、水流缓慢且气候适宜的条件下，藻类易繁殖和停留，湖泊有较高的水华风险，高浓度营养盐也加剧底泥释放的风险。

(a) 2018年1月17日

4 结 论

a.洪泽湖换水周期空间异质性明显，整体呈现出南部湖区水体更新快、北部湖区水体更新慢的特点，换水周期从南到北的梯度变化，主要与洪泽湖的水文情势密切相关。老子山区流场主要受吞吐流影响，换水能力较强，换水周期约为10 d左右；成子湖区流速缓慢，流场主要受风场影响，局部呈现逆向环流，换水能力最弱，换水周期超过100 d。

b. 冬季和春季全湖平均换水周期较长，分别为75 d和60 d，秋季和夏季的全湖平均换水周期较短，分别为49 d和31 d。溧河洼区和临淮区春季较冬季换水能力增强，换水周期分别缩短了31.25%和37.5%。夏季南部湖区的均值小于10 d，受来水量影响，溧河洼区、中心湖区、临淮区和三河区水体更新能力改善显著。

c.实际风场对洪泽湖的换水能力起到正向作用，无风时湖泊水体整体向东南、东北方向流动，换水周期均值为72 d，比无风时减少16 d。风场影响最大的是溧河洼区，比无风时减少89 d，老子山区的水体更新受风场的影响最小。东南风对洪泽湖的换水能力起正向作用，比无风时减少15 d；主要影响溧河洼区，比无风时减少123 d；东北风对洪泽湖的换水能力起到抑制作用，比无风时增加18 d；主要影响临淮区，比无风时增加69 d。