基于二维水流模拟的湖泊型水源保护区划分方法

陶 涛，王韵珏，信昆仑，刘 鹏

（同济大学 环境科学与工程学院，上海 200092）

合理有效地划分饮用水水源保护区，是保护饮用水水源地最大可能免受人类活动影响，保证水质安全的重要措施.国务院2005年印发《关于落实科学发展观加强环境保护的决定》，明确提出“以饮水安全和重点流域治理为重点，加强水污染防治工作.要科学划定和调整饮用水水源保护区，切实加强饮用水水源保护”.因此在传统保护区划分的基础上，研究保护区划分的新方法，为合理保护水源区提供技术支撑尤为重要.本文探讨了基于ArcGIS及二维水流水质模拟的湖泊型饮用水水源保护区的划分方法，以太湖贡湖为研究区域，将望亭立交调水情况作为重要影响因素，应用GIS分析方法和二维水流水质模型，探讨贡湖污染物输移规律，分析湖体水质的变化趋势，研究水源地保护区划分方法，对贡湖饮用水水源保护区做出准确的划分，为贡湖流域的饮用水水源地保护提供切实可用的技术支持.

太湖位于江、浙两省之间，流域面积10000多平方公里，湖面面积2338km2，总蓄水量为44亿m3，年平均吞吐量为52亿m3.近年来因上游补给水不理想，沿湖城市污水排放失控，网箱养殖超负荷以及流域范围内大量使用化肥、农药等原因，近岸水域受到污染.位于太湖流域的无锡市区，自2001年起市区水厂均是以太湖作为供水水源，其中，贡湖水源地承担了无锡市60%的供水，南泉水源厂和锡东水厂的取水点即设置于此，因此贡湖的水质直接关系着无锡市的生活及生产供水能否正常供应.为保障太湖流域供水安全，改善流域水环境，望虞河“引江济太”调水工程于2002年启动.望虞河工程全长60.8 km，包括常熟水利枢纽、望亭水利枢纽及两岸配套建筑物等，工程通过望虞河常熟水利枢纽将长江水引入望虞河，除去望虞河东、西岸引走的水量，其余经望亭水利枢纽进入太湖.“引江济太”调水工程是近期治理太湖的主要应急措施［1］.

1 二维水流水质模拟

1.1 二维水流模型基本方程

水流模型是描述不同水体水文特征和流场时、空分布规律的数学模型.研究不同类型水体流场的时、空分布规律是研究污染物分布的基础，掌握水体中流场特性才有可能研究污染物的分布特性.因此，水流模型是水质模型的流场基础.

本文研究对象贡湖为浅水湖泊，风场、温度等导致垂直混合强烈，因此可假设太湖水体为均匀不可压缩流体，即流体的密度为常值，水体在垂直方向上服从静压力分布.描述湖水运动的二维浅水波控制方程可表示为［2－3］

式中：Z为水位，m；U，V分别为x，y方向上的单宽流量，m3·s－1·m－1；u，v分别为x，y 方向上的流速，m·s－1；g 为重力加速度，m·s－2；h为平均水深，m；c为相应水深下的谢才系数，n 为曼宁系数（河底粗糙系数）；f为科氏系数，f＝2ω·sinψ，ω为地球自转角速度，ψ为湖泊所处纬度；ρ为水的密度；τwx和τwy分别为风应力沿x和y方向的分量，可采用下式计算：

式中：ρα为空气密度；cD为阻力系数；W 为离水面10 m高处的风速矢量；Wx，Wy分别为离水面10m高处的风速沿x和y方向的分量.

1.2 二维水质模型基本方程

描述湖泊污染物输移转化的模型由湖泊水动力模型即水流模型和物质输移对流扩散模型、化学转化模型即水质模型构成.目前在水动力学模拟方面，二维和三维模型的研究都比较广泛，特别是三维模型，已取得了不少研究成果［4－7］.但三维模型比较复杂，成熟应用仍然较少.本研究主要侧重饮用水水源保护区的划分，二维模型分析结果已满足研究需求，因此本文水质模型主要采用二维模型进行分析计算.

在湖泊水质模拟的过程中，本研究主要参照中国水利水电科学研究院建立的污染物质在水体中运动描述的对流扩散的二维方程［8－9］.该模型已成功应用于滇池、太湖、呼伦湖等大型浅水湖泊的水流水质模拟计算，效果较好，本次计算采用该水质模型.二维模型分为平面二维和竖向二维两种情况，前者是指水体的流速和污染物仅在水平面的纵向、横向变化，在竖向（水深方向）混合均匀；后者是指水体的流速和污染物仅在纵向和水深方向变化，横向保持不变.对于浅水湖泊水质问题可简化为水平二维来处理，描述污染物在水体中输移转化的平面二维运动方程如下［10］：

式中：ρ为污染物质量浓度，mg·L－1；Dx，Dy分别为x，y方向上的弥散系数，m2·s－1；μx，μy分别为x，y方向的流速分量，m·s－1；S为各种源汇项之和，g·m－2·s－1；F（ρ）为生化反应项.

参照中国水利水电科学研究院模型，主要水质指标包括高锰酸盐指数（CODMn）、总磷（TP）、总氮（TN）.各污染物指标的生化反应项均作一级简化处理，高锰酸盐指数考虑自净衰减，通过自净衰减系数反映；TP，TN考虑各种因素引起的释放与沉降，通过综合沉降和释放系数反映，不考虑浮游植物对TP，TN的影响，各类水质的生化反应项F（ρ）分别做如下处理：

式中：Kc为 CODMn的自净衰减系数，d－1；ρc为CODMn的质量浓度，mg·L－1；h为平均水深，m；Kp为TP的综合降解系数，d－1；ρp为TP的质量浓度，mg·L－1；Sp为TP的底泥释放系数，g·m－2·d－1；Kn为TN的综合降解系数，d－1；ρn为TN的质量浓度，mg·L－1；Sn为TN的底泥释放系数，g·m－2·d－1.

1.3 二维水流水质模型求解

本研究采用有限元法［11－13］对模型控制方程进行离散处理，主要包括两个部分：控制方程的离散化和代数方程的求解.

网格生成是实现计算区域离散化的关键步骤，采用有限元法模拟实际地形边界，采集研究区域的地形信息资料包括贡湖地理坐标和河底高程，运用TriGrid软件自动生成网格，该软件采用Delaunay三角剖分，将平面区域划分为众多三角单元，见图1.

图1 贡湖水域网格化处理Fig.1 Grid partitioning of Gong Lake

对控制方程离散后得到代数方程，通过高斯消元法即可求得模型数值解，水流模型与水质模型离散格式与网格划分一致，计算同时进行.

高斯消元法求解方程的过程如下.

设有n元线性方程组

高斯消元法由消元和回代两个过程组成，消元的过程就是对方程组（7）的增广矩阵

作有限次的初等变化，使它的系数矩阵变为上三角阵（主对角线以下的元素全部为零的n阶方阵）.经过n－1次消元，得到以 ［A（n），b（n）］作为增广矩阵的上三角线性方程组，见式（10），其与原方程组是同解方程组.回代过程就是由方程组（10）的最后一个方程解出xn，然后通过逐步回代，以此解出xn－1，xn－2，……，x1.

2 望亭立交调水情景分析

根据太湖水体水质特征，选取高锰酸盐指数（CODMn）、总磷（TP）、总氮（TN）作为主要考察指标，水质模拟参数参照中国科学研究院研究成果［14］，见表1.根据太湖水质变化情况和望亭立交调水情况，选择2010年2月和10月作为望亭立交调水的典型月份，因5—9月为太湖流域洪水期，此期间不向太湖调水，故选择5月作为非调水时期的典型月份.利用本文建立的二维水流水质模型，选取3种情景模拟，即2月、10月、5月的数据（见表2）模拟对应日期贡湖流域水流水质变化情况，得出贡湖流场图、流速分布图和水质类别图，以此来分析2010年贡湖水体水流水质变化.

2.12 月望亭立交调水情景分析

以模拟的第5天即3月1日为例，在春季盛行西北风的影响下湖流由贡湖湾顶向贡湖口（东南）方向流动，见图2.由于湖岸的阻滞作用，沿岸湖流流速较慢，同时由于湖中心纵深较大，故此处的流速也有一定的下降，总体呈现如图3的流速分布，湖面平均流速0.023m·s－1.

表1 水质模拟参数Tab.1 Parameters of water quality model

为了防止“大数吃掉小数”［5］，对总磷和总氮指标分别扩大100倍和10倍后加和，绘制2月份贡湖综合水质类别图，见图4.分析得出，对所评价的3项水质指标，高锰酸盐指数以Ⅲ类为主，并在湖中心呈Ⅱ类水质；总磷以Ⅲ类为主，在近岸处降为Ⅳ类水质；总氮以Ⅳ类为主，只在湖中心呈Ⅲ类水质.总体而言，在2010年3月份望虞河水利枢纽调水的情况下，贡湖水质以Ⅲ类为主，其中高锰酸盐指数较好，总氮指数较高.

表2 调水日期及相关数据Tab.2 Dates of water transfer and other data

2.2 10月望亭立交调水情景分析

以模拟的第5天即11月4日为例，在夏秋季盛行偏东风的影响下，湖流由贡湖湾顶向贡湖口（西北）方向流动，见图5.流速分布同样呈现近岸低，由近岸向湖中心流速逐渐增大的规律，总体呈现如图6的流速分布，湖面平均流速0.013m·s－1.

绘制10月份贡湖综合水质类别图，见图7.对所评价的3项水质指标分析得出，高锰酸盐指数以Ⅱ类为主；总磷以Ⅱ类为主，在近岸区域下降为Ⅲ—Ⅳ类；总氮以Ⅲ类为主，在近岸区下降为Ⅳ类.总体而言，对于模拟的在2010年10月望虞河水利枢纽调水的情况下，贡湖水质以Ⅲ类为主，其中总磷和高锰酸盐指数较好.10月贡湖水质明显好于3月的水质，主要原因在于春季是农耕的开始季节，大量化肥使用在春季.据分析太湖湖泊生态系统研究站2008年3月的监测数据，太湖水体溶解性氮平均占水体总氮的79%［15］，且大多数氮以硝酸根的形式存在，因此TN，TP在秋冬季变化趋势明显.

2.3 汛期停止调水情景分析

以模拟的第5天即6月5日为例，2010年5月，当望虞河水利枢纽停止调水时，在夏季盛行偏东南风的影响下，湖流在贡湖内大体上呈现逆时针环流，见图8.流速分布同样呈现了近岸低，由近岸向湖中心流速逐渐增大的规律，总体呈现出如图9的流速分布，湖面平均流速0.014m·s－1.

绘制5月贡湖综合水质类别图，见图10.对所评价的3项水质指标分析得出，5月贡湖水质相对较差，主要原因是望虞河水力枢纽停止调水，贡湖缺少清洁水源的调节，湖面流速较慢，且高温天气适宜蓝藻等藻类的生殖.总体而言，高锰酸盐指数以Ⅲ类为主，总磷以Ⅲ类为主，在近岸区域下降为IV类，总氮以Ⅳ类为主，并在近岸区域下降为Ⅴ类，指标较差.

3 贡湖饮用水水源保护区划分

3.1 一级保护区划分

依据《饮用水水源保护区划分技术规范》（HT／T338—2007），贡湖饮用水水源地一级保护区划分采用经验法，即取水口半径500m内区域划为一级保护区.本文利用上述研究理论，通过模拟2007—2010年水质变化情况，一级保护区内全部水域水质达到Ⅲ级以上标准，其中南泉水源厂一级保护区内Ⅱ类水体占水域总面积不到10%，见图11；锡东水厂在2009和2010年一级保护区内Ⅱ类水体占水域总面积40%以上，见图12.

3.2 二级保护区划分

具体分析贡湖饮用水水源二级保护区，若采用经验法划分，即在锡东水厂和南泉水源厂取水点周围2000m范围定为二级保护区，从模拟的情况上看，该范围内水质均能达到Ⅲ类水体要求（见图13，14），符合饮用水源要求，故采用经验法划分是合理的.

同时分析2010年贡湖流域水质类别图，由于贡湖湖中心水质优于近岸水质，故可对二级保护区进行优化，将保护范围扩大，向湖中心延伸，以提高水源保护区的水质安全性.具体划分方法如下：

（1）利用ArcGIS的等值线生成技术，绘制2010年贡湖流域综合水质指标等值线，根据水质指标等值线划分饮用水水源保护区.

（2）由于南泉水源厂水源地和锡东水厂水源地同属于贡湖流域，距离较近，且考虑到贡湖湖中心水质优于近岸水质，将两水源地二级保护区向湖中心拓展，并合并为一个二级保护区，统一保护.

（3）利用ArcGIS的测距工具，测量岸边污染带的宽度约为2km，考虑到划分的实际可操作性，可将贡湖流域距岸边2km以上的湖面区域全部划定为二级保护区，以提高原水安全可靠性，具体划分结果见图15.

图15 2010年贡湖综合水质指标等值线及二级保护区边界图Fig.15 Contour of comprehensive water quality index and border of the secondary area of water source protection of Gong Lake in 2010

根据模拟结果，优化后贡湖流域内水源地二级保护区内水质基本上呈Ⅲ类水体，Ⅲ类以上水质占保护区水域总面积90%以上，满足饮用水原水水质安全标准，见图16.

图16 优化后贡湖流域水源地二级保护区水质类别Fig.16 Water quality classification in the optimized secondary area of water source protection of Gong Lake

4 结论

（1）系统研究了二维水流水质模型，模拟2010年不同情景模式下贡湖流场和水质指标CODMn，TP，TN的变化过程，模拟结果显示贡湖流场的流速分布总体上呈现出近岸低，由近岸到湖中心流速逐渐增快的趋势，而贡湖水质受流速分布的影响较大，基本上呈现出与流速分布一致的规律，即中心好，近岸水质较差；贡湖水质中，以高锰酸盐指数情况为最好，总磷次之，总氮情况较差.除总氮以外，贡湖湖面大部分区域均能达到Ⅲ类要求.

（2）根据《饮用水水源保护区划分技术规范》（HJ／T338—2007），基于水质模型对水质状况的分析，结合保护区划分原则，并充分考虑2010年内不同调水情景下贡湖流域水流水质变化情况，得到太湖贡湖流域饮用水水源一级、二级保护区范围：一级保护区仍采用经验法划分方法，即南泉水源厂和锡东水厂取水口500m范围内水域划定为一级保护区；二级保护区采用数值模拟法划分，考虑实际操作的可行性，将南泉水源厂和锡东水厂的二级保护区向湖中心延伸并合并，避开岸边污染带，将距离岸边2km以上的贡湖流域划定为水源地二级保护区，保护区面积96.02km2，Ⅲ类以上水体占保护区面积的90%以上.

［1］周小平，翟淑萍，袁粒.2007—2008年引江济太调水对太湖水质改善效果分析［J］.水资源保护，2010，26（1）：40.ZHOU Xiaoping，ZHAI Shuping，YUAN Li.Influences of water diversion from Yangtze River to Tai Lake on water quality of Taihu Lake between 2007 and 2008 ［J］.Water Resources Protection，2010，26（1）：40.

［2］高廷耀，顾国维，周琪.水污染控制工程［M］.3版.北京：高等教育出版社，2007.GAO Tingyao，GU Guowei，ZHOU Qi.Water pollution control engineering［M］.3rd ed.Beijing：Higher Education Press，2007.

［3］CHEN Zhihe，FANG Hongwei，LIU Bin.Numerical simulation of wind－induced motion in suspended sediment transport［J］.Journal of Hydrodynamics，2007，19（6）：698.

［4］LIU Xiaobo，PENG Wenqi，HE Guojian，et al.A coupled model of hydrodynamics and water quality for Yuqiao Reservoir in Haihe River Basin［J］.Journal of Hydrodynamics，2008，20（5）：574.

［5］罗潋葱，秦伯强.基于三维浅水模式的太湖水动力数值试验：盛行风作用下的太湖流场特征［J］.水动力学研究与进展：A辑，2003，18（6）：686.LUO Liancong，QIN Boqiang.Numerical simulation based on a three－dimensional shallow－water hydrodynamic model in Tai Lake—current circulations in Tai Lake with prevailing windforcing［J］.Journal of Hydrodynamics，2003，18（6）：686.

［6］胡维平，秦伯强，濮培民.太湖水动力学三维数值试验研究：3.马山围垦对太湖风生流的影响［J］.湖泊科学，2000，12（4）：335.HU Weiping，QIN Boqiang，PU Peimin.Three－dimensional numerical experiments on hydrodynamics in Tai Lake：3.influence of the reclaimation in Mashan District on wind－driven current［J］.Journal of Lake Sciences，2000，12（4）：335.

［7］逄勇，濮培民.大气－水耦合模式下三维太湖湖流场研究［J］.湖泊科学，1996，8（2）：97.PANG Yong，PU Peimin.Numerical study of wind－driven current in Tai Lake with the air－water coupling model［J］.Journal of Lake Sciences，1996，8（2）：97.

［8］马巍，李锦秀，田向荣，等.滇池水污染治理及防治对策研究［J］.中国水利水电科学研究院学报，2007，5（1）：8.MA Wei，LI Jinxiu，TIAN Xiangrong，et al.Investigation on countermeasures for water environment management and water pollution prevention in Dianchi Lake ［J］.Journal of China Institute of Water Resources and Hydropower Research，2007，5（1）：8.

［9］彭文启，毛战坡，王世岩.洱海、滇池水质保护与水污染防治专题研究［R］.北京：中国水利水电科学研究院，2005.PENG Wenqi，MAO Zhanpo，WANG Shiyan.Water quality protection and water pollution control research in Erhai and Dianchi Lake［R］.Beijing：Institute of Water Conservancy and Hydroelectric Power Research，2005.

［10］宋新山，邓伟.环境数学模型［M］.北京：科学出版社，2004.SONG Xinshan，DENG Wei.Environmental mathematics model［M］.Beijing：Science Press，2004.

［11］周昌玉，贺小华.有限元分析的基本方法及工程应用［M］.北京：化学工业出版社，2006.ZHOU Changyu，HE Xiaohua.The method and application of finite element analysis［M］.Beijing：Chemical Industry Press，2006.

［12］同济大学计算数学研究室.现代数值数学和计算［M］.上海：同济大学出版社，2004.Computational Mathematics Department of Tongji University.Modern numerical calculation and analysis ［M］.Shanghai：Tongji University Press，2004.

［13］YIN Hailong，XU Zuxin，YAO Yijun.Eco－hydraulics techniques for controlling eutrophication of small scenery lakes—a case study of Ludao Lake in Shanghai ［J］. Journal of Hydrodynamics，2007，19（6）：776.

［14］中国水利水电科学研究院.太湖纳污能力分析计算［R］.北京：中国水利水电科学研究院，2002.Institute of Water Conservancy and Hydroelectric Power Research.Calculation of water environment capacity in Tai Lake［M］.Beijing：Institute of Water Conservancy and Hydroelectric Power Research，2002.

［15］朱广伟.太湖水质的时空分异特征及其与水华的关系［J］.长江流域资源与环境，2009，18（5）：439.ZHU Guangwei.Spatio－temporal distribution pattern of water quality in Tai Lake and its relation with cyanobacterial blooms［J］.Resources and Environment in the Yangtze Basin，2009，18（5）：439.