基于MIKE21的南四湖上级湖水量水质响应模拟研究

宫雪亮 ，孙 蓉，芦昌兴，孙秀玲

(山东大学 土建与水利学院，山东 济南 250061)

南水北调工程是一项促进水资源优化配置的重大战略性基础工程。其中，南水北调东线项目成功的关键是污染控制。东线工程第一阶段的主要供水要求满足三类地表水质量要求。南四湖独特的地理位置和东线工程的调水方式决定了当地小流域的水质对保证调水质量起着决定性的作用。

关于水域水质水量模拟的研究起源于20世纪60年代。1965年，Hansne等建立了首个关于浅水海域水位变化的二维函数模型；20世纪70年代，国外学者提出了基于水质水量的湖泊水动力模拟方案，日本、美国等科研机构分别对琵琶湖(日本)、伊利湖(美国)等湖泊水流流态进行了数值模拟，并获得了相应研究成果[1,2]。1975年，Gallagher与Liggerttatal等在前人基础上考虑了气象要素(以风为例)，构建了二维风生流数学模型[3]；1985年，Simon等提出了建立三维模型的开发评价方案；1993年，Kersmiita Zic等改进了该模型，并以日本琵琶湖为例模拟了该地的内波变化以及湖环流结构[4]。2004年，Chao等计算了密西西比流域三维浅水牛扼湖的水动力数值模型，其结果与实际情况吻合度较高[5]。

由于南四湖湖泊属于浅水湖泊，所以本文模拟选取二维模拟，基于丹麦水科所开发的MIKE21软件，该软件主要应用于港口、河流、湖泊、河口及海岸水动力、泥沙及水质的模拟研究[6]。MIKE21中的水动力学模型在计算过程中可以综合考虑各种因素的影响，如地形、气象条件、下垫面情况等，在一定程度上保证了模拟的精确度。

1 构建模型

1.1 研究区概况

南四湖是位于山东省西南部，淮河流域北部的一处狭长水域，处于东经 116°34′～117°21′，北纬 34°27′～35°20′之间，由南阳湖、独山湖、昭阳湖、微山湖这 4 个湖泊相连组成。

南四湖属于浅水大型淡水湖泊，平均水深1.46 m，最大水深2.76 m。1960年在湖腰处建成拦湖大坝，称为二级坝，坝上兴建了节制闸和船闸，二级坝将南四湖分为上级湖和下级湖2个部分。二级坝以北为上级湖，入上级湖河流共29条，集流面积26 934 km2，占总集流面积的88.4%;二级坝以南为下级湖，入下级湖河流共24条，集流面积3 519 km2，仅占总集流面积的11.6%。南四湖流域位于半湿润季风气候区，平均年降水量在750 mm左右，年平均风速为2.7～3.8 m/s。

1.2 确定计算区域

上级湖湖区边界根据经纬度转化为WGS_1984_UTM_ZONE_50N通用投影坐标。根据投影坐标在MIKE软件中得到上级湖的湖区边界，并将此区域划定为计算区域进行网格划分。

本次湖泊模拟所选取的单元网格划分方式为三角形网格，相较于四边形网格划分方式，划分区域更加细致稳定。本次研究上级湖共划分18 188个网格，其中湖中南阳岛为高出水面湖心岛，不进行划分。图1为上级湖计算区域网格详细划分以及局部网格展示。

图1 计算区域网格划分、网格细致图Fig.1 Calculate area meshing, grid detailing

1.3 确定模型边界

根据南四湖上级湖入湖河流，选取汇水面积较大的河流，分为湖东湖西两个区域。湖东区选取6条入湖河流，包括洸府河、泗河、白马河、界河、北沙河、城郭河；湖西区选取7条入湖河流，包括洙水河、洙赵新河、万福河、老万福河、蔡河、惠河、东鱼河。下级湖通过二级坝入上级湖。图2为入南四湖流域划分图。

图2 入南四湖流域划分Fig.2 Watershed division into Nansi Lake

根据入湖河流流域划分情况以及流域汇流面积情况将入湖河流编号整理，形成如表1，最终形成输入模型的边界条件。

1.4 模型初始数据及输入数据

模型初始数据包括处理后的上级湖DEM数字高程数据、上级湖湖区初始水质和湖区初始水位。模型输入项包括各入湖河流的入流流量及污染物数据，出湖水位数据，湖内风速风向以及降水蒸发量。

(1)南四湖上级湖DEM数字高程数据：高程点数据一共包含41 152个原始数据，将高程点数据经纬度转化为WGS_1984_UTM_ZONE_50N通用投影坐标，并对数据中不符合实际情况的个别高程点数据进行相对处理，得到真实可靠的上级湖基本地形数据。

表1 南四湖上级湖入湖河流及编号 km2Tab.1 Incoming rivers and numbers of the upstream section of Nansi Lake

(2)入湖流量数据：假设流域各河流入湖流量总量为100,150和 200 m3/s，由于很多河流未设立水文站，所以使用水文比拟法分配至各入湖河流，进行计算。

(3)风向风速：根据文献资料设置模型上级湖内风速风向，冬季盛行东北风和西北风，夏季盛行东南风和西南风，多年平均风速2.9 m/s，最大风速12.5 m/s。本次研究选取多年平均风速。

(4)湖面蒸发：南四湖湖区面积较大，蒸发量不作为主要变量进行考虑，根据文献资料[7]设置模型上级湖内降水蒸发资料，如表2所示。根据多年平均资料设置不同月份的蒸发数据的时间序列文件。

表2 南四湖多年平均月水面蒸发量表Tab.2 Annual average monthly water surface evaporation of Nansi Lake

(5)本次研究分析湖内水质污染物包括化学需氧量和总磷浓度含量。上级湖湖区初始水质由于监测数据时间尺度较大，加之近几年湖内水质有明显改善，大部分湖区可达三类水水质，所以湖内初始水质定为三类水标准，即COD=20 mg/L，TP=0.05 mg/L。

1.5 模型参数选取

本次湖泊模拟分为水动力模拟与水质模拟两个部分，其中水动力模拟需要控制参数为湖区介质糙率，水质模拟需要控制参数为污染物水平扩散系数以及污染物降解系数。通过对以上参数进行率定，确定最终模型参数。

南四湖湖内情况复杂，航道、深槽、湖草、芦苇、鱼池等遍布湖区，相互交错，且疏密程度有所差异，所以湖区内不同地区的糙率也需要进行概化分区处理。经查阅资料[8]以及分析，最终确定采用4种情况的糙率为：航道(深槽)n=0.025，芦苇n=0.619，湖草n=0.176，明湖n=0.084，,图3为南四湖上级湖湖内粗糙率mesh文件。

图3 湖内介质糙率分布Fig.3 Distribution of medium roughness in the lake

参考文献[9]得出降解系数在0.001～0.1 d-1之间，基于此数值范围参考近年来对于南四湖污染物降解的研究报告，最终确定数值：COD降解系数为0.010 6 d-1，TP降解系数为0.003 1 d-1。

参考文献[10,11]得出扩散系数在0.01～20 m2/s之间，基于此数值范围参考近年来对于南四湖污染物扩散的研究报告，最终确定扩散系数为0.15 m2/s。

1.6 模型验证

1.6.1 水位验证

选取2010年数据进行模型水量水质参数的验证。模型给定计算边界条件及参数选择：梁济运河作为自然入湖河流，上级湖通过二级坝出水渠向下级湖泄水，作为上级湖的出湖口，出湖水位以二级坝实测水位数据为准。

根据上级湖内实测水位资料，选取上级湖内测站南阳站作为模型验证站点，将实测数据与模拟数据进行比较分析。模拟尺度为1 h，输出数据尺度为1 d。

南阳站位于东经116.54°北纬35.01°，处于湖内中游位置，经过对照发现模拟数值与实际水位相差不大，如图4(a)所示，曲线趋势基本相似，波峰与波峰相对，波谷与波谷相对，7-9月份呈现总体上升趋势，10月份之后，非汛期时段内呈现总体下降趋势，与实际情况相符合。南阳站水位相对误差图见图4(b)，整体误差均在可以接受的误差范围之内。实测全年水位平均值为34.305 m，模拟全年水位平均值为34.368 m，相差0.063 m，全年平均水位较为符合。说明建立的模型能够较好地模拟湖区水动力的变化，可以用于各种情景分析。

图4 南阳站水位Fig.4 Water level of Nanyang Station

1.6.2 水质模型验证

水质模型验证是在基于水动力模型模拟基础上，对模型的水质参数进行验证。验证站选取南阳站作为水质模型验证的参证站。选取2010年6月-8月进行模拟，前期为模型预热，最终根据入湖流量以及水质测量时间(7月10日左右)选取7月7日-13日的水质数据与实测数据进行验证。选取COD和TP两种污染物作为水质模拟污染物指标。下面就这两种污染物指标，对模型进行验证分析。

图5为南阳站2010年7月7日-13日的COD和TP浓度模拟值与实测值对比图。污染物浓度总体变化幅度不大，由于汛期有持续入湖径流，所以湖内污染物浓度有上升的趋势，实测值在模拟值变化范围之内。由于总磷易聚集，较难降解的特性，导致在短时间内上升速度较快。综合两种污染物的情况，水质数据均能由模型很好地体现出来，所以认为模型有代表性并且可应用于后续的湖内水质模拟。

图5 南阳站污染物浓度模拟值与实测值对比Fig.5 Comparison between simulated and measured pollutant concentration of Nanyang Station

2 水量水质模拟

2.1 模拟方案设定

二维水质模型是建立在水动力模型的基础上进行模拟分析，本文旨在研究南四湖上级湖的水量水质响应结果，所以方案设置以此为重点。南水北调调水由下级湖向上级湖经二级坝入湖至梁济运河出湖口调水出湖，本文仅对特殊干旱时期调水量在80 m3/s情况进行模拟分析。

湖内初始水位设为兴利水位34.5 m，风向为东北风(315°)，初始水质为三类水。由于二级坝向上级湖调水已经经过下级湖的调蓄过程，所以假设调水水质也已达到三类水，即COD=20 mg/L，TP=0.05 mg/L。

假设流域各河流入湖流量总量为100,150和200 m3/s，使用水文比拟法分配至各入湖河流。入湖径流水质选择较差的Ⅴ类水，即COD=40 mg/L，TP=0.2 mg/L。为了模拟在及时关闸后，湖内总体的污染物变化情况，模型运行时长设定为180 d，其中前20 d入湖河流有入流。

2.2 模拟结果

图6为入湖流量100,150和200 m3/s时梁济运河调水出口断面污染物浓度变化图。可以看出，梁济运河出口断面污染物浓度在开始是持续下降的，这是由于湖内自净使得湖内本身污染物浓度下降；十几天之后会有一个上升，是因为临近梁济运河出口断面的4条河流的污染物逐步扩散至出口断面，使得断面污染物浓度有明显上升趋势；在30天前后达到峰值，之后便持续下降，在50～70 d又开始有小幅度抬升，是由于较远入湖口携带的污染物经过与湖水的充分混合，由湖内水流运动流至出口断面使出口断面有一个较小的峰值，但由于这部分污染物已经经过湖水的充分混合，所以这个峰值会比初始峰值小。

图7为不同情况下上级湖内平均流速以及流向图，在调水情况下，水流流向由南向北，湖内流向为由二级坝调水入湖，以及各入湖口入湖并向湖内扩散，最终经由梁济运河抽水继续北调。流域入湖流量越大，水流运动越快。各方案不同时间水质达到三类水面积占比情况见表3。图8为120 d时，不同入湖流量下湖内污染物浓度情况。

3 结 论

根据模拟结果,可以得到以下结论:

(1)图6、7可以看出,尽管调水出口断面污染物浓度变化过程相似,但流域入湖流量越大，出口断面污染物浓度峰值出现越早,说明在入湖河流来水水质情况一定(Ⅴ类水)的情况下，湖内流速的快慢决定了湖内污染物扩散运移速度。流域入湖流量越大，流速越快，污染物运移速度也就越快，污染物越容易扩散。

图6 不同入湖流量下梁济运河断面污染物浓度变化Fig.6 Pollutant concentration change of Liangji Canal section under different inflow

图7 不同入湖流量下湖内流速分布及流向Fig.7 The velocity distribution and flow direction in the lake under different inflow

表3 各方案不同时间三类水面积占比Tab.3 Water area proportion of third type in different time of each scheme

图8 120 d时不同入湖流量下湖内污染物浓度Fig.8 Pollutant concentration in the lake under different inflow on day 120

(2)虽然入湖流量越大污染物越容易扩散,但来水量越大，也就代表污染物携带量越多,如表3所示,同一时间同一污染物水质达到三类水面积占比,随着入湖流量增大而减小。

(3)当流域入湖流量与入湖水质一定时，不同污染物所表现出的情况也不尽相同。COD降解速率快，最终上级湖内COD达到三类水水质的面积也较多；TP降解速率较慢，水质达到三类水的面积与COD相差较大。

(4)考虑地形因素，由图8可以看出，最终湖内污染物较为聚集的地方均为湖岸边。由于湖岸边水深较浅，水流流速较缓，污染物容易聚集，所以相较于湖内或者是输水航道线上的水来说，最终得到的水质较差。在后续的湖内水资源治理中，可以重点考虑湖岸带的污染治理，会使湖内整体水质上升。

(5)入湖污染物浓度必须严格控制，否则一旦进入到湖内，如果没有大量流域入湖流量或调水流量的冲洗，污染物很难扩散降解，最终导致湖内整体水质下降。尤其是非调水期，由于调水期过后湖内整体初始水位不高，所以需要入上级湖河流开闸放水以补充湖内水量，此时没有调水流量冲洗，就更需要对入湖径流水质进行严格控制，否则入湖径流虽然可以加速湖内污染物扩散运移，但自身所携带的大量污染物质反而对湖内水质改善起到反作用。