引水改善星海湖的水动力-水质特性数值模拟

牛晓瑜, 吴梦迪, 林 玲, 孙秋慧, 徐国宾

(1.天津大学 水利工程仿真与安全国家重点实验室, 天津300354; 2.中水北方勘测设计研究有限责任公司, 天津 300222)

1 研究背景

城市湖泊是城市中重要的生态元素[1],但在城市现代化发展过程中,城市湖泊由于水源单一、流动性差和自净能力低等的特点常常导致其水体富营养化、水质恶化等水环境污染问题。采取水环境改善、引水调度等湖泊综合整治措施对于可持续发展理念下的城市资源规划具有重要作用[2]。星海湖是古黄河游移过程中形成的自然湖泊,经过多年建设,已成为具有防洪和维护贺兰山区域生态平衡[3]等综合功能的城市湖泊。随着星海湖周围的开发利用,入湖污染负荷超过了湖体的自净能力,产生了一定程度的水生态环境污染等问题。星海湖仅由黄河水来补充以提升水质,从而造成了大量的黄河优质水资源耗损。因此,提升湖泊水质与明确引调黄河水量,对保护与节约水资源以及城市的规划发展具有重要意义。

为合理配置水资源以加强保护湖泊湿地等自然资源[4],国内外学者多采用目前较为成熟的MIKE、EFDC(environmental fluid dynamics code)、WASP( water quality analysis simulation program)等模型通过数值模拟方法进行研究[5-11],分析水质恶化的原因,并为实际工程提出相应的治理方案[12-14]。如:Chen等[15]采用EFDC模型研究分析了中国丹江口水库的水质富营养化问题,并对该水库的管理提出了相应的解决措施;李大鸣等[16]建立了MIKE模型对洋河水库进行水质模拟,分析了在水动力作用下污染物随时间迁移扩散的质量浓度变化及分布情况;徐存东等[17]基于MIKE 21模型对沈庄漾河网进行水动力模拟,探明了研究区域水动力提升和流场改善的优化调度方法;鄢碧鹏等[18]采用EFDC模型预测了不同流量、不同出水点与流量分配比例下生态补水对蠡湖透明度的改善效果,为确定蠡湖补水工程建设规模和运行方式提供了依据;王亚宁等[19]以COD和氨氮为输出目标构建太湖EFDC水质模型,探讨了不同入湖水量与污染来源条件下的水质改善率与水质边界敏感性;唐继张等[20]基于MIKE 21建立了二维水动力-染色剂耦合模型,定量模拟分析了西安昆明池(试验段)的换水能力。

本文采用MIKE模型对星海湖进行模拟研究。虽然已有学者[3, 21-22]通过GIS-遥感技术解译影像与水量平衡原理,研究得出了星海湖适宜生态需水总量与在不同水质标准下的生态补水量,并提出了合理的水资源配置方案,但未详细地反映湖泊水动力与水质的时空状况。本研究构建了二维水动力-水质模型,通过模拟计算星海湖不同水文特征、不同引黄河水量的12种工况,分析了各工况的流速、流场、滞水面积、污染物时空分布等水动力与水质特性,明确了最佳引黄河水量,以期为调度黄河水改善星海湖水环境提供理论依据。

2 研究区域概况与研究方法

2.1 研究区域概况

星海湖位于黄河上游宁蒙河段宁夏回族自治区石嘴山市北部贺兰山前洪冲积倾斜平原的下缘,地理位置介于东经105°58′～106°59′,北纬38°22′～39°23′之间,水域面积约为23.42 km2,分为北域、中域、南域、东域和新月海共5个区域,共同承担着大武口沟、归韭沟、大风沟和汝萁沟的来水,星海湖地理位置与湖区分布如图1所示[21]。研究区域属于中温带干旱气候,干旱少雨,日照充足,降雨集中,蒸发强烈,年平均气温为8.7～8.9 ℃,多年平均降水量为179.8 mm,多集中在6—9月,多年平均水面蒸发量为1 291 mm[23]。

图1 星海湖地理位置与湖区分布示意图

星海湖于2004年起通过第二农场渠补充黄河水,年均生态补水量为2 000×104m3。星海湖水质在2020年6月停止补充黄河水后逐渐恶化,7月中域水质为Ⅴ类～劣Ⅴ类[24],水体富营养化。针对星海湖水环境恶化问题,石嘴山市于2021年起对星海湖生态环境进行整治,实施补水、水体内循环、湿地生态修复、污水净化、水资源合理利用、功能区综合治理等多种措施,进一步增强其防洪调蓄和生态安全功能[3]。

2.2 研究方法

2.2.1 MIKE 21 FM模型

(1)水动力基本方程[25]

二维非恒定浅水方程组:

(1)

动量方程:

(2)

(3)

式中:t为时间,s;x,y为笛卡尔坐标系坐标,m;H为总水深,其值为静止水深h与水位ξ之和,m;p、q分别为x、y方向坐标上的流通通量,m3/(s·m);f为Chezy阻力系数,m1/2/s;P为当地大气压强,Pa;Ω为科氏力系数,s-1,Ω=2ωsinφ,ω为地球自转角速度,rad/s,φ为当地纬度;g为重力加速度,m/s2;fw为风阻力系数;ρ为水的密度,kg/m3;W、Wx、Wy分别为风速及其x、y坐标分量,m/s;τxx、τxy、τyy为有效剪切力分量,N/m2。

(2)水质模型方程[11]

h∑Si

(4)

式中:C为污染物浓度,mg/L;u、v为速度分量,m/s;Ex、Ey分别为x、y方向上的扩散系数,m2/s;Si为源汇项,mg/(L·s)。

(3)模块调用

选用MIKE 21 FM中的水动力模块(HD)与对流扩散模块(AD)。水动力模块为MIKE 21的核心基础模块,主要通过对研究区域的地形进行网格化来处理边界条件,选择求解格式、干湿网格条件、河床糙率并考虑降雨蒸发及风场因素、输入源汇作用实现水动力模拟,从而计算得出研究区域的流速和流场。在对流扩散模块(AD)可设定不同类型的扩散系数来反映在不同水动力条件下不同物质的扩散状况。

2.2.2 模型离散求解 模型采用有限体积法将研究区域离散为若干个不规则的三角形和四边形相结合的混合网格,并且保证研究区域的水量与动量守恒。因新月海较小且相对独立,本研究对其不做模拟。对星海湖其他域地形离散共计生成41 258个网格,24 390个节点,星海湖模型计算区域整体网格划分如图2(a)所示,星海湖南域地形较为复杂,局部网格划分如图2(b)所示。

图2 星海湖模型整体网格与南域局部网格划分

2.3 模型边界条件与初始数据

2.3.1 边界条件

(1)水动力模块(HD):确定污水厂、北域(2个)、中域(3个)、南域(3个)进、出水口为开边界,其余为不过流固壁边界,如图3所示。模型计算启用干湿边界,设置干水深hdry=0.005 m,淹没水深hflood=0.05 m,湿水深hwet=0.10 m。

图3 星海湖模型边界条件与监测点设置

(2)源(汇)项:根据星海湖生态环境整治工程和水体内部循环状况,模型中将北域出水点源4引水至东域进水点源1;将东域出水点源2引水至中域进水点源5;将东域出水点源3引水至南域进水点源6;将南域出水点源7引水至中域进水点源8。源汇项与监测点S1(南域)和S2(中域)位置设置如图3所示。

2.3.2 时间步长与初始数据 考虑到模型计算的稳定性和计算效率,库朗数(Courant number)不能大于1;模型计算时长为245 d(2022年3月1日—11月1日),时间步长为600 s,循环计算两次,以保证模型的稳定。

(1)水动力模块(HD)初始数据:本文水文数据均来源于《星海湖生态环境问题整治可行性研究报告》。北域初始水位设置为1 098.33 m,中域为1 098.10 m,南域为1 097.20 m,东域为1 095.35 m,设定研究区域起始为静止状态。星海湖风场设置为恒定值,风速2.45 m/s,风向180°。降水量、水面蒸发量数据采用多年平均数据,详见表1。

表1 星海湖各月多年平均降水量及蒸发量

(2)流扩散模块(AD)初始数据:根据星海湖水质检测结果,选取化学需氧量(COD)为模型水质指标。星海湖北域、中域、南域、东域COD初始浓度分别设定为20、18、18、30 mg/L。

2.4 模型参数率定与验证

2.4.1 参数率定 根据星海湖湖底情况,涡黏系数选用Smagorinsky公式进行估算,取值为0.28。计算床底摩擦力选用曼宁糙率系数n,n取值为0.031 25 s/m1/3。由于水体对污染物有吸附、迁移、沉降与降解等作用,根据资料与模型率定验证,确定星海湖研究区域的COD降解系数为0.002～0.005 d-1。

2.4.2 水动力与对流扩散模型验证 选取2022年5月19日至6月1日星海湖实测数据进行水动力与对流扩散模型验证,验证结果如表2所示。验证结果表明,模型模拟计算出的水位和水质结果与实测值基本接近;星海湖水位模拟计算值与实测值的平均绝对误差为0.02 m,COD浓度模拟计算值与实测值的平均绝对误差为0.69 mg/L。说明本次研究所建立的水动力水质模型精度较高,能够准确模拟星海湖水动力和水质的变化规律,可以用于模拟计算。

表2 星海湖监测点实测数据与模拟计算数据误差统计

3 结果与分析

3.1 模拟计算工况

星海湖需水量整体呈逐年增长势态,为维持周边生态环境,改善星海湖水环境恶化问题,对星海湖进行模拟研究。根据2021年起实施的生态环境整治工程,并结合往年向星海湖引黄河水量的实际情况,在丰水年、平水年、枯水年3种不同水文特征年份下分别设计不引黄河水(引水量为0 )和引黄河水水量分别为700×104、1 400×104、2 100×104m3/a共12种组合工况,如表3所示。

表3 星海湖不同水文特征下引黄河水量模拟计算工况 104 m3/a

工况1、5、9模拟星海湖未实施黄河补水工程的水动力水质情况;工况2、6、10模拟向南域进水口引黄河水100×104m3/a,向中域两个进水口各引黄河水300×104m3/a,再由北域出水口排出700×104m3/a水量;工况3、7、11模拟向南域进水口引黄河水200×104m3/a,向中域两个进水口各引黄河水600×104m3/a,再由北域出水口排出1 400×104m3/a水量;工况4、8、12模拟向南域进水口引黄河水300×104m3/a,向中域两个进水口各引黄河水900×104m3/a,再由北域出水口排出2 100×104m3/a水量。

3.2 模拟计算结果分析

3.2.1 水动力分析 各个工况下星海湖水动力模拟计算结果如表4所示。

表4 不同工况星海湖水动力模拟计算结果

由表4可知:(1)星海湖引水稳定后,12种工况中流速波动幅度较小,南域和中域平均流速在0.006～0.008 m/s之间,南域、中域最大流速分别在0.065～0.072 、0.114～0.118 m/s之间。(2)在相同水文特征年份中,丰水年工况9、10、11、12的南域滞水区面积比例较枯水年工况1、2、3、4相应削减了5.89、3.90、3.75、3.80个百分点,较平水年工况5、6、7、8相应削减了0.92、1.21、0.78、1.07个百分点。在相同引水量工况下,不引黄河水工况1、5、9的南域滞水区面积比例较引黄河水700×104m3/a工况2、6、10相应增大了3.39、1.11、1.40个百分点,较引黄河水1 400×104m3/a工况3、7、11相应增大了4.33、2.33、2.19个百分点,较引黄河水2 100×104m3/a工况4、8、12相应增大了5.11、2.87、3.02个百分点。(3)各个工况中域滞水面积比例在11.70%～12.20%之间,相差较小。中域滞水面积比例比南域小36.49～44.12个百分点。上述结果表明,不同水文特征对星海湖南域滞水面积比例的影响大于引不同黄河水量对其的影响,但两者对中域滞水区面积比例的影响均较小。星海湖水域面积大,湖内整体流速较低,不同工况下星海湖整体及星海湖中域局部流场均呈现相同趋势,故以工况3(枯水年引黄河水1 400×104m3/a)为例进行流场模拟,结果见图4。

图4 星海湖整体流场及中域局部流场模拟结果(工况3)

由图4(a)可知,星海湖水体由南向北流动,由于湖中岛屿众多,在引水稳定后出现了复杂的水流流态,最大流速仅出现在入水口与边壁位置。图4(b)显示,星海湖中域北部沿边壁形成了较大的逆时针环流,流速呈由内向外增大的趋势,环流中心流速基本低于0.002 m/s,沿岸流速最大为0.021 m/s。中域南部远离引水主路线区域流速较低,同时因地形影响出现多个逆时针小环流,建议在类似区域增设曝气机以增大流速,改善局部水动力,提高水体循环能力。

3.2.2 水质分析 以星海湖实测结果为模型的初始条件,模拟得到不同水文特征年份与不同引黄河水量工况下星海湖水质指标COD浓度随引水历时的变化趋势,如图5所示。由图5可以看出,不同模拟工况中星海湖水质指标COD浓度变化趋势相同,虽然水质改善效果略有差别,但均能够满足对星海湖水质在Ⅳ类及以上的要求[26]。星海湖南域COD浓度在0～50 d内先小幅度下降后再升高,50～170 d内再急速下降后又缓慢升高,170 d后缓慢下降,总体而言处于小幅波动变化状态。对比不同工况下星海湖南域水质随引水历时的变化,工况2、3在引水170 d后,COD浓度达到最高,分别为20.52、20.04 mg/L,分别有87.35%,96.56%的时段能够达到地表水Ⅲ类水质标准。其他工况下水质均能达到地表水Ⅲ类水质标准,工况3、6、9分别在枯水年、平水年、丰水年中为最优方案。建议根据COD浓度变化规律分时段引黄河水。星海湖中域COD浓度在0～50 d内小幅度波动,50 d后急速下降再缓慢升高,总体变幅较小。由于中域来水量大,能够提高水体流动速度,从而加快了污染物COD的降解。星海湖中域在12种工况下水质均能达到地表水Ⅲ类水质标准20 mg/L。

图5 星海湖南域和中域COD浓度随引水历时变化趋势

星海湖引水稳定后南域和中域的COD浓度及其削减率计算结果如表5所示。由表5可知,星海湖南域COD削减率在枯水年工况1～4和平水年工况5～6为负,分析其原因,虽然引黄河水能够有效降低星海湖COD浓度,但因其南域整体流速低,地形变化大,且枯水年与平水年来水较少,导致该区域滞水面积比例较大,水体降解能力不足。整体结果表明,星海湖南域丰水年各工况COD平均浓度比枯水年、平水年各工况COD平均浓度分别降低了12.73%、7.36%,相应的COD削减率平均分别提高了13.79和7.63个百分点。星海湖中域丰水年各工况COD平均浓度比枯水年、平水年各工况COD平均浓度分别降低了4.95%、2.72%,相应的COD削减率平均分别提高了4.80和2.64个百分点。在同一水文特征年各工况中,随着引水量的增加,南域和中域的COD浓度逐渐降低,削减率逐渐提高;在相同引水量工况中,丰水年的COD削减率普遍高于枯水年与平水年,丰水年引水对星海湖水质的改善效果优于枯水年与平水年。

表5 不同工况引水稳定后星海湖南域和中域的COD浓度及其削减率

引水稳定后,在12种工况下星海湖水质空间分布情况相似,故以工况2为例进行分析。不同引水天数星海湖COD浓度空间分布如图6所示。

图6 不同引水天数星海湖COD浓度空间分布(工况2)

由图6可知,由于大风沟、大武口沟和污水厂等来水导致污染物COD主要聚集在星海湖南域和北域的西北部。南域岛屿众多,水体流速低,污染物容易积累于岛屿周围,COD平均浓度为19.87 mg/L;中域因地形与流场因素,环流中心处COD浓度较高,平均为16.64 mg/L,边界COD浓度较低,平均为13.58 mg/L;北域水道崎岖且流速低,导致COD容易积累,建议增设人工湿地以降解污染物。

4 讨 论

本研究所构建的星海湖水动力水质耦合模型能够较为全面地反映星海湖水动力特性与水质指标COD的时空分布,确定了在不同水文特征年星海湖最优的引黄河水量。李添雨等[14]同样采用此类模型量化分析了沙河水库水量水质变化情况,也取得了良好的模拟效果。

4.1 水动力时空变化特征分析

本研究结果表明,由于引水口分布的影响,星海湖在入水口、边壁位置与靠近引水主路线区域流速较高,远离引水主路线区域流速较低,环流流速由外向内逐渐降低,这会造成星海湖水动力不均衡,水体得不到充分流动,难以达到全湖水体流通的效果。这与李悦等[13]对大通湖采用MIKE 21模型进行研究时得出的结论相同。此外星海湖南域因岛屿众多,水体流通不畅,流速较低,导致滞水面积增大,建议设置多个补水口或增设曝气机以改善星海湖水动力情况。

4.2 水质时空变化特征分析

本研究结果表明,星海湖水质总体呈现上升—下降—上升—缓慢下降的变化趋势。由于星海湖承接汝箕沟、大风沟、归韭沟、大武口沟等沟道洪水,春季冰雪消融导致来水较多,并且水中大量藻类及附近水产养殖生物开始大量繁殖,导致水体污染加重。夏季降雨集中且洪水同步形成,同时气温升高导致蒸发量增大进一步使得水环境恶化,污染物浓度升高。本文研究结果与王世强等[21]对星海湖需要加大补水的月份的结论相吻合,但本次模拟结果能够更加直观地反映污染物的时空分布情况。因此建议在春、夏季加大对星海湖的水环境治理,合理调控黄河水资源的分配。

5 结 论

本研究建立了星海湖水动力水质耦合模型,模拟计算了星海湖在不同水文特征年、不同引黄河水量的12种工况,以流速、滞水区面积比例和COD浓度削减率等作为评价指标对星海湖水动力与水质改善效果进行了分析,得到以下结论:

(1)星海湖水动力与水质改善效果与不同水文特征及不同引黄河水量相关:在相同水文特征工况中,增加引黄河水量能够进一步减少滞水面积比例,增加污染物削减率,水质指标COD浓度逐渐降低,星海湖水动力与水质改善明显;在相同引黄河水量工况中,由于丰水年来水量多于枯水年与平水年,星海湖水动力状况改善明显,从而使得丰水年水质指标COD削减率与滞水面积削减比例普遍高于枯水年与平水年。

(2)考虑经济性与水质改善效果,建议在枯水年、平水年、丰水年分别选择工况3、6、9引水方案。由于春季和夏秋季来水量大,COD浓度增高,建议考虑分时段引黄河水。

(3)星海湖水域面积大,流场复杂且流速较低,滞水面积占比较高,水体难以自循环,可考虑局部增设进水口与曝气机,加强水体流动性。星海湖水体中COD主要来源于污水厂、大风沟、大武口沟等来水,建议在入湖口采取污染物削减措施,以降低其入湖浓度,并考虑增设生态岛、人工湿地等设施增强水体降解能力。