基于MIKE21数值模拟的滆湖换水能力评估分析研究

Evaluation and Analysis of Water Exchange Capacity of Gehu Lake Based on MIKE21 Numerical Simulation

摘要: 湖泊换水能力直接影响着污染物迁移和水体自净能力，为探究不同引水流量入湖对滆湖换水能力的影响，文中使用MIKE21建立了滆湖水动力水质模型，模拟并对比分析多种工况及工况组合对湖区整体和局部的水体交换能力的影响；提出了使用水体更新时间和水体交换率两个指标，定量分析交叉验证水体交换能力。结果表明：引水流量为100、150、200、250 m³/s时整个湖泊水体交换率到达80%以上分别需要60 d、45 d、32 d、30 d。本文的网格划分方案为湖泊水动力水质模型建模提供了新思路，换水能力两个指标能够从时间和空间上综合全面的评估湖泊换水能力，研究成果可为滆湖水质改善工程提供重要科学参考。

Abstract: The water exchange capacity of lakes directly affects the migration of pollutants and the self-purification capacity of water bodies. To explore the influence of different water diversion flows on the water exchange capacity of the lake, MIKE21 was used to establish a hydrodynamic water quality model of the lake, and the influence of various combinations of working conditions on the overall and local water exchange ca-pacity of the lake area was simulated and compared. The quantitative analysis of water renewal time and water exchange rate was proposed to cross-verify the water exchange capacity. Considering the invest-ment amount and effect, the results show that when the water diversion flow is 100, 150, 200 and 250 m³/s, the water exchange rate reaches over 80% of the whole lake needs 60 d, 45 d, 32 d and 30 d, respec-tively. The meshing scheme provides a new idea for the modeling of lake hydrodynamic water quality model, and the dual indicators of water exchange capacity can comprehensively evaluate the lake water change capacity both from time and space, which can provide an important scientific reference for the wa-ter quality improvement project of Gehu Lake.

1. 引言

湖泊在水文循环过程中扮演着重要的角色，是地球上大气圈、岩石圈、生物圈、水圈等各个圈层相互作用的重要节点。近几十年来城市建设进程迅猛，很多城市湖泊受到建设活动和周边用地的影响，出现了面积萎缩、水动力条件不足、水质恶化等一系列问题 [1] 。现有研究通常使用引水入湖方式来提升水体水动力条件，增强湖泊自我修复能力，改善湖泊水质水环境状况 [2] 。

换水能力并没有一种统一的指标来定义，一般的评价指标有水体交换率、水体更新时间、水龄 [3] 、滞留时间、半交换时间 [4] 、净化时间、水体存留时间 [5] 等。黄春琳等 [6] 使用水龄分布特征来表征太湖水体交换能力，黄爱平等 [7] 建立鄱阳湖水龄模型描述其换水能力，唐继张 [8] 等以换水周期为指标观察昆明湖换水能力。当前关于湖泊换水能力的研究通常都是单一指标分析，或从空间上分析湖泊各局部区域的换水特性，或从时间上分析湖泊整体浓度变化，难以兼顾湖泊换水能力在时间和空间、局部和整体的变化特点。

为探究引水入湖对滆湖水动力条件的影响，定量分析引水流量变化对滆湖换水能力的改善程度，本文提出水体更新时间和水体交换率两个指标来评估湖泊换水能力，既能兼顾其时空变化和局部与整体之间的关系，又能使两个指标的结果形成交叉验证提高结果的可信度。水体更新时间是指水体内物质量减少到原有物质总量的37%时所需的时间，常用于评估近岸海域水交换的研究，近几年被广泛应用到湖泊、河流等水体的运动特性和物质输移的研究，能够表现出水体各个区域水体更新的空间分布特征 [9] 。水体交换率可表征整个水体中某种物质含量的变化，体现整个水体中物质平均浓度随时间的变化情况。使用水体更新时间和水体交换率两个指标来描述滆湖这样较为封闭湖泊的换水能力，同时研究其局部水体更新特征和整体水体交换能力是更加全面的。

目前湖泊水动力模型包括MIKE21、EFDC、Delft3D、太湖模型 [10] [11] [12] 等国内外知名的水质水动力模拟软件。本文采用MIKE21软件，建立滆湖水动力水质模型来研究其换水能力。采用对流扩散方法，以水体更新时间和水体交换率两个指标来定量评价换水能力，探讨引水入湖情况下滆湖的水动力条件和换水能力的变化。研究结果对水质改善工程规划具有的实际参考价值，为湖泊水体更新研究提供理论依据。

2. 研究区域概况

滆湖位于江苏省常州市西南方向，横跨常州市武进区和宜兴市，属于浅水湖泊，平均水深1.3 m，最大水深4 m，湖泊长度22 km，最大宽度9 km，现状水域总面积约为146 km²，总库容2.1亿m³ [13] [14] 。如图1湖底地形3D图所示，湖泊形状西岸岸线较为平滑，东岸岸线较为曲折，相比湖泊西岸，其死水区域较多。湖泊上设有5个观测点(观测点位置见图1)观测1#断面处引水流量为100 m³/s，2#、3#断面处水位为0.2 m时，5个观测点坐标及各点处稳定流速度和水深见表1所示。此次模拟设置1个规划进水口(见图1，断面1#)和2个规划出水口(见图1，断面2#、3#)。

Figure 1. 3D map of bottom terrain in the Gehu Lake

图1. 滆湖湖底地形3D图

Table 1. Coordinates, stable flow velocity and water depth of five observation points in the Gehu Lake

表1. 滆湖5个观测点坐标和稳定流速及水深

3. 研究方法

3.1. 水动力模型构建

1) 数据预处理与模型建立

滆湖水动力学模型采用丹麦水利研究所(DHI)开发的MIKE21数学模型软件进行建模，MIKE21模型是二维平面水流数学模型的一种。该模型在国内已经得到了广泛的运用，常用于河流、水库、湖泊、海湾等水体区域的流场、水质和泥沙运输的模拟。其详细的原理与数学求解方法在本领域很多研究中都有极为详尽的介绍 [15] [16] ，本文就不再过多介绍。

数据预处理过程中，利用CAD和ArcGIS软件将测绘得到的湖泊边界及高程散点数据的dxf文件转换成xyz文件并导入MIKE21中生成模型。将湖泊5个观测点的经纬度坐标转换到相应的投影坐标下，利用MIKE21自带的数据转换工具进行批量转换并生成观测点坐标文件导入模型。

2) 网格划分及网格划分方案对比

为适应湖泊复杂的湖底地形和曲折的岸线情况，本文在建模过程中绘制了4种不同特点的网格划分方案(如图2、图3所示)，来对比不同网格划分方法对模型模拟精度和计算效率的影响，并比选出最合适的网格划分方案。考虑到湖泊引水出入口处流态更加复杂，所以在B、C、D方案中在出入口处不同程度的增加了计算网格的

Figure 2. Mesh diagram of A and B mesh division schemes

图2. A、B网格划分方案网格图

Figure 3. Mesh diagram of C, D meshing schemes

图3. C、D网格划分方案网格图

Table 2. Attribute tables for four meshing schemes

表2. 四种网格划分方案的属性表

密度；同时考虑到引水出入口处流向更加确定，所以C方案中在左出口(2#断面)处设置了矩形网格，D方案中出入口(1#、2#、3#断面)处均设置了矩形网格。4个方案网格划分主要特点如下：A、B方案均采用三角网格划分，A方案网格划分在全湖区域内较为均匀，而B方案在入水口和出水口处网格划分进行了有明显的细化处理；C、D方案采用湖泊中心使用三角网格、出湖口的河道区域使用矩形网格，D方案相比C方案的网格更加密集。各方案的网格形态、网格参数设置见表2。

3) 模型主要参数设置

模型重要参数及模型边界设置信息见表3。

Table 3. Construction and parameter setting of MIKE21 hydrodynamic model

表3. MIKE21水动力模型的构建和参数设定

3.2. 模型率定与验证

1) 最优网格方案选取

分别对4种网格方案进行多次率定，得到每个网格划分方案下的最优糙率空间分布图，将相应的最优糙率分布图导入对应网格方案模拟计算，得到各方案观测点流速和水深的模拟值，进行误差分析和比较。由分析表4和表5可知，A方案5个观测点水深平均相对误差为3.4%，相对误差最大值为7.00%；B方案水深平均相对误差为0.44%，相对误差最大值为0.75%；C方案水深平均相对误差为3.70%，最大相对误差为9.4%；D方案水深平均相对误差为1.75%，最大相对误差为3.90%。4个方案水深模拟结果最大相对误差均在10%以内，5个观测点平均误差均在4%以内，模拟结果均在合理范围之内。

Table 4. Relative errors of water depth simulation values for each scheme

表4. 各方案水深模拟值相对误差

而流速的相对误差各方案之间差距较大，A、B方案平均误差分别为25.66%、90.76%，误差过大，与实测值差距过大；D方案平均误差16.75%，但最大误差为47.19%，与实测值仍有较大差距；C方案平均误差12.26%，最大误差不超过25.0%，流速模拟精度与实测值相差不大，可认为C方案模拟得到湖泊流场与实际情况基本相符。

Table 5. Relative error of flow velocity simulation value of each scheme

表5. 各方案流速模拟值相对误差

在相同计算机配置、相同计算时长和步长情况下，各方案的运行时长见表6所示。从表6可见，在相同计算机配置和相同计算时长、计算步长情况下C方案的运行时间最短，说明C方案的计算效率最高，同时也说明计算网格数量越少计算效率越高。

综合各网格划分方案的水深、流速模拟精度、计算效率，可认为四个方案中C方案为最优的模型网格划分方案，并使用C方案作为最终模型方案。

Table 6. The running time of each program

表6. 各方案的运行时长

2) 最优方案的结果

通过不断地对模拟期水动力参数进行率定和调整，得到C方案的最优糙率分布图如图4所示。湖泊糙率分布设置靠近湖中心位置糙率较小，岸边区域靠近死水区、淤塞区域糙率较大。进行水动力模拟得到滆湖流场矢量图如图5所示，模拟结果与实际流场基本一致。

从图6、图7可见，模型验证结果表明，经过一定时间的模拟模型的水深、流速均收敛到稳定值，并在小幅度范围内波动与实际情况相符。在模型运行过程中，水深在模型模拟9天便达到稳定，各观测点水深分别收敛至1.115 m、1.089 m、1.326 m、0.885 m、2.546 m，各观测点水深与对应实测稳定水深平均相对误差在3.70%以内。流速在运行过程中波动较大，在图6中已剔除前25天流速波动过大的数据。流速在45天左右达到稳定，各观测点流速分别收敛到0.096 cm/s、0.042 cm/s、0.096 cm/s、0.045 cm/s、0.757 cm/s，各观测点流速与对应实测稳定流速平均相对误差在12.26%以内。

模型验证结果表明，该模型能够较为准确的模拟滆湖水动力情况，可以用于水体更新时间和水体交换率的研究。

Figure 4. Optimal roughness distribution map of scheme C Figure 5. Illustration of flow field vector in the Gehu Lake

图4. 方案C最优糙率分布图 图5. 滆湖流场矢量图

Figure 6. Variation curve of water depth at each observation points Figure 7. Variation curve of flow velocity at each observation points

图6. 各观测点水深变化曲线 图7. 各观测点流速变化曲线

3.3. 换水能力两个指标定义

1) 水体更新时间

水体更新时间是基于水箱模型的一种水交换定义指标 [5] ，是用于评估湖泊水体内纯粹由于对流扩散作用导致湖泊物质交换速率的一个常用指标。可以理解为湖泊内均匀分布某种稳定不可降解的溶解态示踪剂。本文设置这种示踪剂浓度为1 (无量纲)，从湖泊入口处引入清水且示踪剂浓度为0。湖泊水体引水过程中，受对流扩散作用影响水体中各网格内示踪剂浓度降低到初始浓度的37%时(即示踪剂浓度降低至0.37时)，所消耗的时间为水体更新时间。水体更新时间表达式如下：

$T_{renewal}=T\left(C\left(t\right)\le C_0\times 37\%\right)$ (1)

式中： $T_{renewal}$ 表示水体更新时间， $C\left(t\right)$ 表示t时刻示踪剂浓度， $C_0$ 表示示踪剂初始浓度， $T\left(C\left(t\right)\le C_0\times 37\%\right)$ 表示示踪剂浓度降低到初始浓度37%以下时的时间。

2) 水体交换率

水体交换率是指湖泊水体内部浓度为1的与上述相同的稳定示踪剂，在经过湖泊对流扩散作用一段时间后，输移至水体外部的物质总量与初始物质总量的比值 [17] 。水体交换率反映湖泊内水体整体的平均换水程度。水体交换率的表达式如下：

$E=\frac{C_0\cdot S_0\cdot h-{{\displaystyle \sum }}_n^N{C}_n\left(t\right)\cdot S\left(n\right)\cdot h}{C_0\cdot S_0\cdot h}$ (2)

式中： $E$ 表示水体交换率， $C_0$ 、 $S_0$ 、 $h$ 分别表示初始浓度、湖泊总面积、湖泊平均水深， $N$ 表示湖泊网格总数为1461个， $n$ 表示网格序号， $C_n\left(t\right)$ 第 $n$ 个网格 时刻的示踪剂浓度， $S\left(n\right)$ 第 $n$ 个网格的面积。

3.4. 模拟方案与模型工况设置

1) 模拟方案

本文以水体更新时间和水体交换率作为指标来评估研究滆湖水体交换能力。模拟方案为：如图1所示的1#断面作为湖泊引水入口，2#、3#断面为出水口，在引水入流的水力驱动下改善湖泊水动力条件，增加湖泊水体对流扩散能力。为方便进行水体更新时间和水体交换率的计算，模型假设湖泊水体充分混合，并存在某种不可降解的稳定溶解态示踪剂浓度为1 (无量纲)，引水清水入流示踪剂浓度为0。

2) 工况设置

湖泊当前原工况引水流量为100 m³/s，相较于引水前湖泊的水动力条件和水质已经得到了较为明显的改善。但由于湖泊死水区域较多，特别是由于右岸区域岸线曲折导致右岸区域水体更新依然较慢。因此，为进一步提高湖泊水动力条件、对比和量化不同引水流量对湖泊水体交换能力的影响，使模拟结果能更好的服务于实际引水工程改造，本文设置了3种加大引水流量的工况。工况1：加大引水流量至150 m³/s；工况2：加大引水流量至200 m³/s；工况3：加大引水流量至250 m³/s。

4. 结果与分析

4.1. 水动力特征

引水流量为100 m³/s时，根据图5流场矢量图显示，滆湖河口处流速较大且流速方向较为一致，在右侧出口处岸线凹槽处流向较为紊乱。在湖体中央流速相较湖边流速更大且方向更加有序，说明引水入湖能够有效的提高湖泊的流动性，增强水动力条件。两岸流速相较于湖中心流速较慢，且流向相对更加无序，右岸部分区域形成涡流，说明虽然湖泊的水动力条件增强但对湖泊的死水区域影响较小，仍需进一步加大流量或采取其他措施来增强右岸区域的水体流动性。流场整体状况与实际情况相近，说明模拟结果具有可靠性。

4.2. 水体更新时间

根据滆湖水动力模型模拟结果绘制的各工况下水体更新时间图9可见，湖泊的水体更新时间大致存在以下规律：越靠近出入水口，水体更新时间越短，反之则越长，且湖泊东岸存在较多死水区；湖泊中央区域较左右岸边区域的水体更新时间更短。

根据图8可见，引水流量加大后湖泊中心区域、左岸的水体更新时间均缩短，右岸水体更新时间过长的区域面积有明显变小，水体更新时间小于10天的区域明显变大。根据UNEP流域综合管理指南建议湖泊水体更新时间一般不超过30天。当入水口流量达到200 m³/s时，湖泊由中心向两岸扩散的区域整体水体更新时间基本下降至30天以内。而引水流量由200 m³/s增加至250 m³/s时，工况3相比于工况2水体更新时间在30天以内的区域面积并没有很明显的增加，在右岸区域水体更新时间大于60天的区域虽然稍有减少，但依然存在死水区域

Figure 8. Water renewal time chart of water body under various working conditions

图8. 各工况水体更新时间图

水体更新时间过长的问题。加大入流流量能够缩短该湖泊的水体更新时间，但该湖泊右岸部分死水区仍然存在换水周期较长，污染物滞留的现象。

4.3. 水体交换率

根据滆湖水动力水质模型各工况下的示踪剂变化模拟结果，计算得到滆湖15 d、30 d、45 d、60 d、75 d、90 d的水体交换率。滆湖水体交换率计算结果表见表7，水体交换率变化曲线见图9。

Table 7. Water exchange rate under various working conditions

表7. 各工况下水体交换率

Figure 9. Variation diagram of water exchange rate under various working conditions

图9. 各工况水体交换率变化图

由图9可见，各工况的水体交换率随引水时间的延长而增大，且均呈现前期水体交换率的增长较快，后期水体交换率的增长较慢的趋势。说明引水时间越长湖泊内原有的污染物因对流扩散作用运移到湖泊之外的污染物越多，且越到后期因湖泊内污染物浓度的降低水体交换率也随之变化得越慢。将四个工况作比较发现，引水流量越大水体交换率的增长越快，相同时间内水体交换率变化也越大，说明加大引水流量是加快降低湖泊整体污染物浓度的一个有效办法。

从表7中展示的数据可以看到，经过15天的引水入湖，原工况和工况1、2、3的水体交换率分别达到了24.67%、36.86%、49.52%和60.01%。说明在短时间内，加大引水流量能够快速明显提高水体交换率，当引水流量达到200 m³/s以上时，滆湖水体能够在15天就能将湖泊污染物整体浓度降低到原来一半的效果。进一步分析可以看到，水体交换率达到80%以上，原工况需要60天、工况1需要45天、工况2需要32天、工况3需要30天，说明引水流量越大整体浓度降低80%以上所需的时间越短，与水体更新时间计算结果形成了相互验证。

从湖泊水动力条件来看，湖泊流场与实际情况相符合，湖中心流速较大两岸流速较小，水动力特征明显。从水体更新时间来看，工况2引水流量为200 m³/s时湖泊中心区域大部分水体更新时间均已经达到30天以内，湖泊右岸局部区域水体更新时间过长；且工况3继续增大引水流量未见特别明显改善。从水体交换率来看，湖泊整体的水体交换率达到80%以上工况2需要30天左右，工况3需要30天。

5. 结论

1) 在建模划分网格时，加密出入水口处网格的密度能够降低水深模拟的相对误差，但对流速模拟的相对误差没有明显提升作用，甚至可能会增大误差。在二维湖泊水动力学模型当中，河道区域使用矩形网格或者三角网格对模型模拟精度没有明显影响。模型网格的数量能够明显影响到模型运行速度，所以在保证精度的情况下可以尽量减少网格数量以提高模型运行效率。

2) 模型模拟结果水深和流速均能够收敛到稳定值，且模拟值与实测值相对误差较小、结果相吻合，说明所构建的模型具有合理性。

3) 增大湖泊引水流量能够明显提高湖泊的水动力条件。当工况2流量增加到200 m³/s时湖泊大部分区域水体更新时间均小于30天，湖泊水体15天交换率可达到近50%、30天湖泊水体交换率达到近80%。因考虑到工况3继续加大补水流量带来的边际效益递减，所以本文建议在实际工程中，引水流量加大至200 m³/s即可。既可以缩短湖泊各区域水体更新时间，又能提高湖泊整体的水体交换率，保持湖泊水质指标的稳定。

4) 同时使用水体更新时间和水体交换率两个指标来评估湖泊的换水能力是更加系统和全面的。水体更新时间可以更加清晰直观的分析湖泊换水能力空间上的分布特征，能看到湖泊内局部区域乃至每个网格的水体更新时间。水体交换率可帮助分析湖泊换水能力的时间变化特点，能够分析每天乃至每个小时湖泊整体的污染物减少量。两个指标结合，既能看到湖泊换水能力的空间特征，又能掌握湖泊换水能力的时间变化特点，同时观测局部变化又能把握整体情况。

基金项目

国家自然科学基金区域创新发展联合基金(U20A20317)。