基于MIKE21的西安昆明池(试验段)换水能力特征研究

唐继张， 周维博， 安宝军， 杨 浩， 李 颖

(1.长安大学 水利与环境学院， 陕西 西安 710064； 2.长安大学 旱区地下水文与生态效应教育部重点实验室，陕西 西安 710064； 3.陕西省西咸新区沣东新城斗门水库建设管理中心， 陕西 西安 710086)

1 研究背景

城市化进程的快速推进使得现代社会对自然的改造力度越来越大，城市建设发展过程中对水资源的要求也在不断提高[1-4]。近年来，人工湖泊在城市的规划定位中突显了其重要价值[5-8]：不仅具有旅游价值，推动当地及周边的商业发展、提高居民的生活质量，而且一些湖泊还兼具调洪蓄水、保障城市供水安全、改善局地小气候等生态价值。但与天然湖泊相比，人工湖泊较为封闭，水环境单一，且水体流动性较差，纳污及自净能力不足[9-10]，加之周边人类活动的影响，容易引起湖泊富营养化、水质恶化等问题。若不及时有效地进行预防治理，会对城市的建设发展起到阻碍作用。明确人工湖泊的水流场特征及水体交换能力，对预防湖泊水质恶化与污染治理、促进城市水环境健康发展意义重大[11-13]。

许多学者针对湖泊水体交换能力问题进行了大量研究。例如，许莉萍等[14]运用数值方法建立水动力和对流扩散模型，研究了人工湖的水体更新时间和水体交换率。Delhezé等[15]建立了一维河流和二维河口的耦合模型对斯凯尔特河河口对流扩散进行研究，得到了河口水体滞留、输运及影响时间。王冼民等[16]利用EcoTaihu模型，在不同蓄水量方案下研究太湖的氮、磷浓度最低值的换水周期。李云良等[17]运用数值模拟方法，定量对比了不同季节时鄱阳湖的换水周期，得出了鄱阳湖在四季的换水能力大小。Umgiesser等[18]建立了地中海泻湖三维数学模型，得到了潮流流场、水体交换率、水体输运时间和水体混合等情况。李大鸣等[19]基于物质迁移和对流扩散理论，研究了洋河水库中TN、TP、NH3—N质量浓度的变化和分布。姜恒志等[20]针对普兰店湾海域在潮流作用下的水交换能力，将海湾分为不同区域建立潮流数值模型，得到各区域的水体交换率。鉴于数值模拟方法在湖泊研究中的成熟运用[21-22]，本文采用MIKE软件，在水动力基础上耦合染色剂模型，从时间及空间尺度上，对昆明池(试验段)进行模拟研究，结合湖区流场特征，定量分析整个湖区的换水能力，探究湖区在正常运行条件下最可能出现水质恶化问题的水域及原因。研究结果为昆明池(试验段)水体的污染预防与治理提供了科学指导，也为后期昆明池扩建工程提供了理论依据，还对于一般人工湖泊的水环境综合管理及调度运行具有重要理论意义和参考价值。

2 资料来源与研究方法

2.1 研究区概况

西安昆明池(另称斗门水库)位于东经108°45′～108°48′，北纬34°11′～34°13′之间。地处关中平原中部，秦岭以北，属暖温带半湿润大陆性季风气候区，其利用原昆明池遗址进行扩建，分为南湖和北湖，是一座以调蓄“引汉济渭”工程来水、兼顾城市生活供水、沣河防洪及改善生态环境的综合性平原水库工程。

本次研究对象为昆明池试验段工程(以下简称为昆明池)，湖区设有5座人工岛屿，水域面积0.47 km2，库底高程396.8 m，正常蓄水位400.47 m，蓄水深度3.67 m，蓄水量为155×104m3。湖区位于沣河与太平河之间，由沣惠渠引水，再退水至太平河，水系较为丰富。但沣河及沣惠渠水质较差，且围绕昆明池打造的休闲公园和周边建筑群以及一系列的人为活动使得湖区现状水质较差，局部水域出现富营养化状态。研究概况见图1。

图1 研究区概况

2.2 数据来源

本文所用昆明池地形及水位数据来自《陕西省斗门水库工程环境影响报告书》及《西咸新区沣东新城斗门调节水库岩土工程勘察报告》；降雨蒸发、风速风向资料来自于课题组于研究区设置的小型气象站自动监测采集。

2.3 研究方法

2.3.1 MIKE21 FM水动力模型

(1)模型基本方程。MIKE21 FM为平面二维水流数学模型，该模型遵循Navier-Stokes方程，且服从Boussinesq假定和流体静定假设。其二维水流连续方程为：

(1)

二维水流动量方程为：

(2)

(3)

(4)

(5)

考虑流体黏性影响，水流的涡黏项用Tij表示，其值按下式估算：

(6)

(7)

(8)

式中：A为黏性系数，m2/s。

(2)方程组离散求解。模型采用有限体积法进行研究区域的离散求解，离散为若干个不规则的三角形网格，可采用三角形与四边形相结合的混合网格对局部区域进行加密，以便更好地反映复杂边界，有限体积法保证了水量和动量在计算域内的守恒。本次对昆明池地形离散共计生成4 178个网格，2 307个节点，最终剖分网格见图2(图中已标注进水口、出水口及龙头监测点位置)。

(3)模型模块的选择。水动力模块(HD)是模型模拟的基础，主要通过对研究区域的地形网格化处理、控制网格干湿条件、选取合适的求解方法并考虑模型的输入输出作用，从而计算出模拟的研究区水深和流场分布规律。对流扩散(AD)模块通过设置不同类型的扩散系数模拟物质在水体中经过对流和扩散过程而发生的扩散现象特征，其优点在于忽略了污染物质在水中运移时发生的复杂的物理、化学、生物及生态过程，突出了区域水体流场的作用。考虑湖区的水体环境及研究目的，本次模拟选用水动力模块和对流扩散模块作为研究手段。

图2 研究区网格剖分

(4)定解条件的确定。湖岸及湖心岛边界面的法相流速为零，设为零法向流速的陆地边界；进水口为流量，设定为实测的随时间变化的流量边界。出水口为水位，设定为实测的水位边界；模型输入项主要还包括湖区风速风向及降雨蒸发，皆设定为在空间上恒定，随时间变化。

2.3.2 水体交换研究方法

(1)换水周期。关于湖泊水体交换能力的讨论，许多学者从不同角度进行了定量的研究，提出了水龄、水力滞留时间、更新时间、曝光时间和换水周期等概念。根据研究对象及模型特点，本文选用换水周期来定量研究湖区的换水能力，其中换水周期为可溶性保守物质初始的单位浓度经过水流对流扩散作用，衰减到某一程度所需要的时间，表征湖泊对于湖区水体自我净化能力的大小[17,23]。基于湖区内保守物质的最终浓度几乎不可能为零这一事实，本文换水周期选用浓度变化的指数衰减函数表示，具体公式如下：

Ct=C0·e-t/Tf

(9)

式中:t为时间，d;C0为保守物质初始浓度值,g/L;Ct为保守物质在t时刻的浓度值,g/L。

由公式(9)知，当t=Tf(V/Q)时，浓度衰减至初始浓度的e-1，即将剩余浓度降低至初始浓度的37%时所需要的时间定义为换水周期。

(2)频率分布曲线。频率分布曲线[24]能够直观反映整个湖泊的换水周期面积分布，基于频率曲线的定义，通过分析变异性，进一步探讨湖泊的换水能力，评价换水周期的优劣。曲线越陡峭，说明湖泊换水能力变异性越弱，系统换水越快；曲线越平缓，则说明湖泊换水能力变异性越强，系统换水越慢。

2.3.3 模拟方案及模型设置 本文选定MIKE21 FM的水动力和对流扩散模块，在湖泊水动力模型基础上通过耦合染色剂模型，模拟计算出昆明池的流场特征和换水周期分布。根据上文方法定义，本文将昆明池水体设置为单位浓度为1的保守型示踪剂作为初始浓度场，湖区进水口设为浓度为0的流量边界，湖口开边界设为自由出流边界，随着模型的运行，通过监测研究区内每一个网格单元的剩余浓度变化，得到昆明池的换水周期空间分布。

考虑到湖区在夏季降水及蒸发量较大，且周边公园日开放时常较长，人流较多，这一时期湖水水质遭受外界因素的影响较大，故结合现有数据资料，设定本次模拟时段为夏季，即2018年6月1日至9月1日。

2.3.4 模型率定与验证 根据实测水位数据对模型进行反复率定，调参得到的模型主要参数为：干湿网格临界水深为0.05 m、Smagorinsky系数为0.28、曼宁糙率系数为28 m1/3/s。再利用研究区中龙头监测点7月份实测水位数据对模型进行验证，验证结果如图3所示，模拟过程水位与实测基本一致。

图3 龙头监测点水位验证

(10)

依据公式(10)计算得到该模型的Skill值为0.94，大于标准值0.65，表明该模型拟合精度极好，可用于对昆明池的模拟计算。

3 结果与分析

3.1 水动力场模拟

运行模型对昆明池进行水动力模拟，得到夏季昆明池流场模拟结果，在模拟期内，各时段流场变化较为一致，故选择模拟期中间时段某一天(2018年7月26日)作为代表性流场进行分析，结果见图4。由图4可见，昆明池在人工岛附近均形成多个较为明显的环流，方向呈逆时针流动，其成因或为湖区在东侧进水口引入水流，西南侧出水口排水，人工岛在水流运动路径上与其发生碰撞，改变了水流方向形成环流。在沿岸也出现许多较小环流，这可能是由于昆明池湖岸线起伏，特别是北岸在码头、观水台等曲折较大的临水地带，水流受其扰动作用，形成范围较小的环流。由于整个湖区水域在南北方向上呈“中间窄、两边大”的特点，水流在东西方向交换过程中，经过湖中间地带时受到地形阻碍，改变了水流方向，故在湖东西部区域分别形成较大的环流。

图4 2018-07-26昆明池流场模拟结果

受环流等影响，昆明池流速也具有明显分区域的特征，流速分布模拟结果见图5。

由图5可看出，贴近人工岛的环流流速较大，约为0.04～0.06 m/s，但湖东北部人工岛附近，因为水域偏离水流主方向，故流速较低，约为0.002～0.005 m/s。沿湖岸线受地形因素形成的各小环流的近岸流速较大，约为0.04～0.07 m/s，个别区域流速达到0.09 m/s。而在湖中心区域，水流较为缓慢，流速基本低于0.004 m/s，且距离湖岸较远区域流速约小于0.002 m/s。

3.2 换水周期模拟

对昆明池进行换水周期模拟，得到夏季昆明池各区域的换水天数。其中，龙头、出水口处的浓度-时间曲线见图6，分析图6得出，龙头处的换水周期为27 d，出口处换水周期为36 d，即使新水至东侧进水口引入，也分别需要约27、36 d才能对龙头和出水口处水体进行较完全的置换。这可能是由于龙头位置处于地形夹角处，水流循环受阻，而出口位置距进水口较远，水流路径较长，故更新缓慢，换水周期较长。各“时间-浓度”为指数衰减形式，进一步表明了本文所用的换水周期计算方法在本次模拟中具有良好的适用性，体现了模拟结果的可靠、合理性。

图5 2018-07-26昆明池流速分布模拟结果

图6 昆明池龙头和出口站点浓度随时间变化曲线

进一步对整个湖区换水周期分析(见图7)，得出湖区换水周期具有明显的空间异质性，根据其换水周期在东西方向上阶梯递增的特征，大致将湖区划分成4个区域(Ⅰ、Ⅱ、Ⅲ、Ⅳ)。其中：靠近进水口的区域Ⅰ，其水流路径最短，故换水周期最小，约为1～10 d，但东北角区域换水周期有所增长，约为10～17 d，或是因为该区域偏离进水口主流方向，且水流与该处的人工岛发生碰撞，阻碍了其向东北角扩散，故降低了东北角水流交换能力，造成换水周期延长；位于龙头前端的区域Ⅱ，其换水周期约为17～28 d，且与区域Ⅲ在湖岸地形收缩处形成一明显分界线。结合湖区流场特征，分析可能是由于区域Ⅱ存在的较大环流造成的，且环流边界恰位于分界线附近，故在此处，部分水流前进方向改变被迫回流至区域Ⅱ内，形成一明显的换水周期分界线，且随着水流路径的增长，换水周期较区域Ⅰ也有所增大；穿过东部湖区，水流到达西部湖区；区域Ⅲ换水周期约为28～32 d，在小环流，特别是靠近湖南岸一带的小环流作用下，水流流动方向改变，做循环流动，换水周期普遍较大，并且在西侧与区域Ⅳ形成一明显界限；区域Ⅳ处于昆明池最西侧，其进水水流路径最长，普遍换水周期约为32～35 d，在人工岛附近及靠近湖岸一带，部分区域换水周期超过35 d，这同样是受人工岛与湖岸形成的环流影响所致。

图7 昆明池换水周期空间分布

3.3 频率分布曲线

根据昆明池换水周期模拟结果，做出其换水周期频率分布曲线，见图8。分析图8可知：曲线上存在3个明显拐点(分别对应换水周期为10、17、34 d)，将曲线分成4段，呈现出“平缓-陡峭-平缓-陡峭”的特点，表明昆明池换水能力变异性不具备明显的强弱特性，即昆明池不能简单地定义为一个快速换水或慢速换水的湖泊系统。其中，约10%的昆明池湖区换水周期≤10 d，约40%的昆明池湖区换水周期≤17 d，约70%的昆明池换水周期≤34 d，还有不足10%的昆明池湖区换水周期超过了35 d。

图8 昆明池换水能力变异性分析频率分布曲线

整体上看，昆明池流场及换水周期分布均主要受湖中人工岛及湖岸地形因素的影响。在水动力交换弱和受地形影响形成环流的区域，水流更新时间较长，换水周期较大，这种换水周期分布特征与水动力场的密切相关性，表明了模型模拟结果的合理性。通过现场调研，昆明池在人工岛附近及湖西部区域，更容易出现水体富营养化等问题，这也与昆明池模拟结果相印证，进一步说明了本次模型模拟的科学性与结果的合理性。通过对昆明池流场特征及换水周期分布的探究，能够明确昆明池湖区范围内易受水体交换更新不利影响而发生水质恶化问题的区域，有利于对湖区水质问题作出针对性的预防治理措施，维持湖区水环境的健康可持续运行。

4 结 论

(1)昆明池流场受人工岛及湖岸地形影响较大。在各人工岛附近，均形成多个逆时针流向的小环流，近岛处流速较大，平均可达0.04～0.06 m/s；靠近湖岸，存在许多明显小环流，流速约为0.04～0.07 m/s，最大流速可达0.09 m/s；在湖东西部均存在较大环流，流速基本低于0.004 m/s。

(2)昆明池换水周期具有明显空间异质性，在东西方向呈阶梯递增：区域Ⅰ约1～17 d，区域Ⅱ约17～28 d，区域Ⅲ约28～32 d，区域Ⅳ约32～35 d，个别区域超过35 d。

(3)昆明池属于快速换水与慢速换水并存的湖泊系统，约10%的水域换水周期≤10 d， 40%的水域换水周期≤17 d，70%的水域换水周期≤34 d，不足10%的水域换水周期超过35 d。

(4)加强对湖区西部的监测，提高该水域的水体交换能力，是进一步改善昆明池水质，预防湖区出现富营养化问题的关键。