水动力过程对湖泊中物质与能量的迁移转化起到决定性作用, 是湖泊中污染物、泥沙等输移的重要基础; 湖泊水动力特性主要与湖泊本身轮廓特征、吞吐流、风驱动力等物理参数有关(李一平等, 2012, 2014; Yao et al, 2019; Chen et al, 2020)。对于大型浅水湖泊而言, 风和吞吐流为其主要的驱动力, 风场会改变湖水的运动轨迹, 决定着湖泊的环流形态、结构及湖流流速的大小, 影响物质在水体中的迁移与扩散, 进而对湖泊水环境产生影响(姚静等, 2016; Liu et al, 2019)。吞吐流也是湖水运动的主要驱动力之一, 它通过水体交换, 进而对湖泊营养盐扩散与迁移、泥沙冲淤及湖水污染物浓度场的变化产生影响(Gao et al, 2015; 赵正文, 2018)。

在大型浅水湖泊水动力研究过程中, 对于湖泊的实时监测不易实现, 故数值模拟是研究湖泊水动力的重要手段之一。已有研究对湖泊水动力的数值模拟开展得比较深入: 如针对太原市迎泽湖的富营养化问题, 武春芳等(2014)通过调整补水方式、补水频率、湖泊柔性结构三方面, 提出改变水动力条件的方案并进行数值模拟, 结果表明: 通过加大湖泊的水力循环能力、改善入湖口水质的控制方案可以有效改善水质; 鄱阳湖存在倒灌现象, 唐昌新等(2015)借助EFDC的染色剂模块及水龄模块模拟了鄱阳湖的全年水动力变化过程, 长江倒灌的条件是临界来流大于实际总来流, 2008年鄱阳湖的倒灌入流扩散面积占全湖总面积的0.97%, 6.6%, 3.5%, 12.6%, 鄱阳湖水龄可由长江的顶托作用时间加倒灌水体泻出时间估算; 北方湖泊冬季出现结冰现象, 陈友媛等(2012)针对这一现象在东昌湖水动力模型中引入了冰盖因素, 分析六种风应力拖拽系数对水动力的影响, 讨论了冰盖面积对水动力的影响, 发现覆盖冰盖的水域, 风应力对水动力作用会被减弱甚至抵消。但上述研究通常将风场、吞吐流等物理参数作为次要影响因素, 对于风场、吞吐流与大型浅水湖泊流场的响应关系鲜有研究。

乌梁素海是我国典型的大型浅水湖泊, 其南北长35—40 km, 东西宽5—10 km。湖区面积多年维持在293 km, 库容在2.5亿—3亿m³之间。湖区总体水深浅, 平均水深在0.8—1.5 m之间, 最大水深只有3 m左右。湖底高程变化不明显, 南北方向呈现北高南低的趋势, 淤泥层厚, 深水区多为主航道(李兴, 2009; 田伟东等, 2016)。乌梁素海地处西风地带, 风期较长, 风速较大, 年平均风速在2.5—3.4 m/s之间, 年最大风速18—40 m/s。同时乌梁素海也是河套灌区重要的组成部分, 是确保黄河内蒙古段枯水期不断流的重要调蓄补给库, 承接河套灌区百分之九十以上的农田灌溉排水、生活污水、工业废水, 具有改善水质, 减少污染物质对黄河水质直接影响的作用(赵锁志, 2013)。但由于水体流速缓慢、换水周期长、自净能力差, 导致乌梁素海水环境受到了一定的影响(田伟东, 2016)。想要研究乌梁素海的水环境问题, 需要找到根本原因及关键影响因素。水动力特征与规律的研究对改善乌梁素海水环境问题具有重要意义, 精准地模拟和评估风场及吞吐流在湖区的影响形式, 对于湖泊水动力的研究尤为关键, 因此进一步开展乌梁素海水动力的模拟具有宝贵的实际应用价值和科学价值。本文基于二维水动力模型, 设置定常风场和不同算例的吞吐流量, 模拟不同情况下的乌梁素海流场形态, 分析风场及吞吐流对乌梁素海的主要影响区域、影响量级及流态响应机制, 为研究乌梁素海污染物的迁移转化等研究提供理论支撑。

1 研究方法及方案设计 1.1 模型控制方程

水动力模型遵循Navier-Stokes方程、Boussinesq假设和流体的静定假设。基本方程式如下:

式中: t为时间, 单位为d; x, y为笛卡尔坐标系下坐标, 单位为m; η为表面水位, 单位为m; d为静水水深; h为总水深, h=η+d, 单位为m; f是科氏力系数, f=2ΩsinΦ, Ω为地球自转角度速率, Φ为当地维度, 单位为1/s; g为重力加速度, 单位为m/s²; ρ₀为(淡)水的参考密度, 单位为g/ml; ρ为水的密度, 单位为g/mL; τ_sx, τ_sy, τ_bx, τ_by分别为x和y方向的表面风应力和底部切应力, 单位为N; S为源项, 单位为kg/(m³·s); v_s、u_s为源项水流流速, 单位为m/s; u, v为流速在水深上的平均值, 单位为m/s; T_xx, T_xy, T_yy为侧向应力, 单位为N; p_a为当地大气压强, 单位为Pa; s_xx, s_xy, s_yx, s_yy为辐射张量, 单位为Pa。

1.2 水动力模型建立

乌梁素海水动力模型的构建首先需要制作网格文件, 插入湖底地形高程, 然后将网格文件导入水动力模块, 同时输入边界条件、初始条件、相关参数等, 最后通过水动力模块输出结果得到模拟数据(如水位、水深、流速等)。数据的收集整理、模型的构建均由作者完成。

1.2.1 计算区域和网格划分

乌梁素海主要补给源为总排干渠、八排干渠、九排干渠、十排干渠, 其中总排干渠约占总补给量的88%, 湖水由乌毛计退水口排出。进出口位置、湖底高程如图 1所示。乌梁素海平均水深0.8—1.5 m, 相较于平面尺度很小, 故假设水体垂向掺混均匀, 选用垂向均匀的平面二维不可压缩模型。乌梁素海每年11月初至翌年3月末为冰封期, 由于冰盖作用, 风场对流场影响无需考虑, 故本模型的模拟时间为2018年4月1日—11月1日, 时间步长为3600 s。针对乌梁素海湖岸沿线蜿蜒曲折的特点, 选择三角形非结构网格, 三角网格划分区域更加细致稳定, 能够更好地拟合地形和湖岸边界。模型共布置网格单元12005个, 节点数为6643个, 网格尺寸最小为150 m, 最大为280 m, 平均尺寸为240 m。生成网格后, 将乌梁素海湖底的坐标及高程值制作成X, Y, Z格式的文件(X为横坐标, Y为纵坐标, Z为高程值), 输入到模拟软件中, 生成水下地形图。

1.2.2 定解条件

本次模拟的乌梁素海水动力模型上游入湖源汇项包括总排干渠、八排干渠、九排干渠、十排干渠, 下游出湖源汇项为乌梁素海乌毛计退水渠, 各源汇项流量数据分别为各排干渠及乌毛计退水渠2018年4月至10月逐日流量数据。根据《乌梁素海综合治理规划》, 到2020年, 乌梁素海水位控制在1018.8—1019.2 m, 故模型运行的初始水位取1019.2 m(水位值和湖底地形高程值都基于海平面)。乌梁素海位于干旱半干旱地区, 太阳辐射强、蒸发量大、降雨量小、干湿期差异大, 精准的蒸发量和降雨量数据使模型更加稳定, 降雨量和蒸发量数据采用气象站2018年4月至10月每15 min一次的监测数据。风场数据来源与时间区间同上。

1.2.3 参数设置

本模型主要率定的参数为糙率系数, 糙率系数是水力学计算的关键参数, 与水深、床面形态、植被条件等因素有关(吴乔枫等, 2018)。以芦苇为主的多种挺水植物大量生长的区域面积和湖泊明水面积相间分布是乌梁素海主要特征之一, 水面有一半被芦苇覆盖, 所以糙率系数并不是单一的数值。本模拟将遥感影像图在Arcgis中划分为不同的区域类型, 勾勒出湖区现状的各种区域的轮廓图, 对芦苇区、明水区、水道等分别定义, 设置不同的糙率系数(赵振兴, 2013)。相对于比插值法, 这种方法更加准确地定义了湖区不同区域的糙率系数, 更适用于乌梁素海的水文特征, 使模型更加精确。糙率系数值如表 1所示。

1.2.4 模型的验证

由于缺乏乌梁素海湖区实测数据, 故采取人工监测的手段。在湖区设置了12个监测点位(Y1—Y12), 监测频率为每月3次, 监测内容为湖区12个点位的水深, 水深值采用SR-100型湖泊水深测深仪测量, 点位分布如图 1所示。为了减小湖泊水深值测定的结果因环境、仪器操作原因产生的误差, 每个点位测量3次, 取平均值作为该点的实测水深值。本研究采用2018年4—10月湖内12个样点的实测数据和模型运行结果数据进行对比验证, 结果表明: 各点位水深相对误差均在±3%之内, 确定性系数和纳什效率系数也均达到0.96以上。依据点位区域分布特点选取4个代表性点位, 分别为西大滩区Y2、东大滩区Y4、芦苇区Y7、南部明水区Y11, 制作水深验证图(图 2)。

1.3 风场分析及模拟方案

风场数据采用位于坝头村的乌梁素海湖泊湿地观测研究站实测的2015—2018年逐日每15 min一次的风速、风向资料。乌梁素海湖区面积较大, 湖区不同位置风场存在差异, 但由于缺少多个监测站的实测风场资料, 故本研究假定全湖风场在空间分布上一致。本文主要研究乌梁素海盛行风场对湖泊的水动力场影响, 故对于收集到的风场资料进行统计, 得到出现频率较高的风向并同时计算出其对应的平均风速, 模拟该风向、风速对湖泊水动力场的影响方式和差异。通过监测站数据风玫瑰图(图 3)可知, 乌梁素海4—10月以SSW向风场为主, ENE向风场次之, SSW向风场平均风速为2.46 m/s, ENE向风场平均风速为2.56 m/s, 故分别模拟2018年4月1日—11月1日期间无风、SSW向定常风、ENE向定常风3种情况下的水动力变化过程, 2种风向平均风速均取2.5 m/s。风场是影响湖泊水动力的主要因素, 湖泊水动力影响湖泊污染物的迁移、水质的分布等。设置定常风场主要针对真实风场中出现的频率较高的、作用时间较长的风场, 模拟得到乌梁素海湖区大部分时间流场的形态进行研究, 分析水动力对乌梁素海水质分布的影响, 为乌梁素海的水质研究提供一定的理论支撑。

1.4 水情分析及算例设计

根据2018年4—10月各排干渠数据(图 4)得到此期间进出口最大吞吐流量和最小吞吐流量。由于九、十排干渠吞吐流对湖区南部区域影响较小, 故主要研究吞吐流对总排干渠和八排干渠口附近区域的影响。为了对比不同吞吐流量作用下乌梁素海的水流特征, 研究吞吐流与风场2种驱动力的作用关系, 设计了以下4种算例进行模拟, 具体算例如下:

算例1: 总排干渠入湖流量为47.81 m³/s、八排干渠入湖流量为4.4 m³/s, 乌毛计退水渠出湖流量为49.82 m³/s, 无风;

算例2: 总排干渠入湖流量为16.6 m³/s、八排干渠入湖流量为0.7 m³/s, 乌毛计退水渠出湖流量为19.31 m³/s, 无风;

算例3: 总排干渠入湖流量为16.6 m³/s、八排干渠入湖流量为0.7 m³/s, 乌毛计退水渠出湖流量为19.31 m³/s, 2.5 m/s的SSW向定常风场;

算例4: 总排干渠入湖流量为47.81 m³/s、八排干渠入湖流量为4.4 m³/s, 乌毛计退水渠出湖流量为49.82 m³/s, 2.5 m/s的SSW向定常风场。

2 结果与分析 2.1 风场对湖区流速的影响

为了分析风场对乌梁素海流速的影响, 将有风作用下的湖泊流场与无风条件的流速分布进行比较。在无风状态下, 湖区流速分布不均匀, 西大滩区流速较大, 而南部湖区及北部小海子人工湿地水流几乎静止, 造成这种现象的主要原因是: 乌梁素海湖底地形整体地势较为平坦(北高南低), 没有大的湖底地形高程差, 不会因为高程差产生较大的流速; 总排干渠、八排干渠位于湖区的西部, 排干渠的水流进入湖体, 水量越大, 流速越大; 湖泊中部是芦苇茂盛地区, 芦苇面积占25%以上, 水流经过芦苇区, 流速衰减, 随水流方向到达南部时, 流速衰减殆尽。北部小海子同样是芦苇茂密区, 并且北部地势略高于南部, 流速衰减更加明显。由此可以看出, 单一的吞吐流无法满足乌梁素海湖区良好的水动力条件, 风驱动力是乌梁素海流场的关键驱动力, 风场作用下, 大部分湖区流速增加, 湖区明水面出现多处环流, 芦苇区的流速受风场影响较小。

2.2 定常风对湖区流场的影响

乌梁素海芦苇区面积较大, 湖泊被芦苇区分成了几处明水区域。明水面由于没有芦苇遮挡, 出现了受风场影响明显的环流区域。鉴于风场对湖区不同空间的影响形式和效果存在差异, 划分出湖区风生流较为明显的3处区域进行研究对比。图 5中A区域为西大滩区, B区域为东大滩区, C区域为南部湖区。定常风场作用下, 流场于2018年4月1日6:00便趋于稳定状态, 流场形态基本不发生变化。八月相较于其他月份为丰水月, 模型中环流形态效果呈现更加明显, 故选择2018年8月1日10:00的瞬时流场作为典型流场进行分析。

A区域(图 6)位于湖区西部, 来自总排干渠、八排干渠的水流汇入该区域, 在无风条件下, 该区域流速受排干渠吞吐流影响相对较大, 水由排干渠口进入该区域导致进口附近流速较大, 流速随着水流方向递减。在SSW风向作用下, 总排干渠进水口浅水区流向东南方向的水流受到影响, 改变了原本的流向, 向东北方向偏移, 与该区域远离排干渠口的中心区域向西南方向的补偿流形成小尺度顺时针环流, 环流中心流速近乎为0 m/s。A区域东南方向靠近芦苇带的浅水区的水流也受到影响, 流向与风向一致, 流速较大, 超过0.01 m/s, 与中部区域的补偿流及排干渠口附近的水流同时作用形成不规则的"几"字型。在ENE向风场下, 入湖口附近的浅水区水流仍然受到风场的影响, 整体流向东南方向。东南部靠近芦苇带的浅水区域水流流向西南方向, 拉动中部形成补偿流, 区域整体水流呈"S"型, 区域北部伴有小尺度逆时针环流。2种风向在A区域产生的环流形态呈中心对称。

水流从A区域流向B区域(图 7)主要通过连接2个区域的航道实现, 航道内部流速相较于芦苇区较大, 在无风条件下, 由于芦苇的阻挡, 到达B区域后流速逐渐衰减, 水流到达B区域时东部时, 流速衰减为0 m/s, 出现大面积滞流区。在SSW向风场作用下, B区域东南方向浅水区域受到风向的影响, 沿着B区域东南部边缘, 顺着风向流向北部, 拉动中部深水区域形成逆时针环流; 西北部浅水区域水流也受到风场影响, 水流向北方, 与中心补偿流同时作用, 形成2个顺时针的环流; 相较于SSW向风场, ENE向风场下B区域流场变化更加复杂, 由于受到A区域东部水流的影响, 连接A与B区域航道水流流速要高于SSW向风场时的流速, B区域的西北部流场受到航道内水流方向的影响, 大环流被分割, 形成了3个相近的小环流, 东南部由于水深较浅, 加之岛屿掩蔽作用影响, 形成顺时针环流, 环流速度大于该区域其他环流。对比2种风向对该区域的影响, 风场对流场的作用形式相似, 方向相反。

C区域(图 8)是乌梁素海较大的明水面区域, 地形呈中间低四周高的趋势, 无风条件下, 水流经过湖区中部芦苇区后, 流速极小, 近乎停滞。在SSW风场作用下, 东南部及西北部湖滨区沿湖岸方向流动, 由于空旷水面无芦苇阻碍, 沿岸流速较大, 经过中部芦苇区的水流顺着C区域北部的航道汇入中部补偿流中, 形成了2个方向相反的大尺度环流, 西南部近岸区域水流受地形影响, 形成了小尺度逆时针环流。ENE风向对该区域影响形式与SSW风相类似, 同样是东南部及西北部湖滨区流向变化带动中心形成环流, 但由于东方向风向占主导, 导致形成的环流形状不规则, 面积较大, 流速较小。X点处由于芦苇遮挡, 几乎不受风场影响, 但是水流方向依然发生偏转, 导致西北部形成的环流区域远小于东南部环流, 这说明湖泊的边界形态也是改变水流方向的重要因素之一。

2.3 吞吐流对湖区流场的影响

吞吐流是湖水运动的主要驱动力之一, 是影响湖泊水动力的关键因素, 吞吐流量的变化对湖区流场的影响也比较大(姜恒志等, 2013; 吴青, 2017)。在得到各算例的流速矢量图(图 9)后, 对比算例1和算例2, 湖区在仅增大吞吐流量情况下, 湖区的水流形态基本未发生改变, 只是水流的速度发生变化, 吞吐流量增大, 湖区流速增加。对比分析算例2和算例3, 在湖区常规的吞吐流作用下, 加上定常风场驱动, 湖区的水流结构发生明显变化, 这说明湖区的水流形态受风场的影响较大, 单一常规的吞吐流量不会使湖区形成环流; 对比算例3和算例4, 发现增加了吞吐流量后, 湖区整体流速增大, 除了A区域靠近总排干渠口处的环流形状受到轻微影响, 大部分环流区域位置及环流形状不发生改变, 但环流的强度随着吞吐流量的增加而减弱。这说明在吞吐流量增大过程中, 吞吐流驱动力制约了风生流驱动力, 导致环流强度减弱。

3 讨论

本文利用相关资料和数据建立了乌梁素海水动力模型, 针对乌梁素海湖区出现的高频率风向, 模拟得到乌梁素海定常风作用下湖区的流速变化、环流的特点、形态等, 同时研究了不同算例下流场形态与流速的变化情况。定常风场对于湖区流场作用的形式是单一的, 然而真实的风向、风速是随时间空间持续变化的, 并非如假设情况是一成不变的, 故实际的流场也是复杂多变的。通过实际风场作用与假设的定常风场作用下的流场形态进行对比可以得出: 随着时间推移、风场及吞吐流的变化, 实际流场的运动主要影响因子是风场, 吞吐流对湖泊的流场影响相对较小。本研究涉及的SSW、ENE定常风场下的流场形态只是实际情况中出现的风向中的2种, 乌梁素海在不同风场(实际风场)转换过程中, 没有因为风场的转换而影响流场产生特殊的变化(即不同的风场在稳定状态下会对应特定的流场形态)。实际流场的研究基本可以认为是不同时刻的风场对湖泊流场的影响研究。出现这种现象的主要原因是乌梁素海湖底地形平坦, 吞吐流量变化小, 即便是丰水期, 吞吐流量也不会有较大的变化, 不会影响湖区环流的位置, 形态等。已有研究表明(韩龙喜等, 2011), 湖泊风生流过程分为3个阶段, 为湖水流向与风向基本一致阶段、流向逐渐偏离风向的过渡阶段、形成环流的稳定阶段, 风场的持续时间也会对湖泊流场造成影响, 但是本文侧重研究特定风向、吞吐流量对湖流瞬时流态的影响, 极端情况下的风速和吞吐流量对流场的影响以及风场持续时间对流场的的影响本研究没有涉及到。

以往类似风场和吞吐流对湖泊流场影响研究也有很多, 如韩龙喜(2011)对高盐度艾比湖进行三维的风流数值模拟, 结果表明艾比湖沿岸水深较浅的水体受风应力的影响较大, 湖水流向与风向基本相同, 然而水深较大的区域容易形成与风向相反的补偿流。姚静(2016)模拟了无风、NE和SSW 3种定常风向情景下的鄱阳湖水动力模型, 结果显示同一风场作用下, 不同湖区的环流结构、影响形式、持续时间也不尽相同, NE向风场和SSW向风场产生的环流位置相近, 尺寸相近, 但是方向相反。在对狭长型雁鸣湖水动力模拟过程中, 赵正文等(2018)设置了几种吞吐流量来研究湖区流场的变化, 发现随着流量的增加, 湖区环流位置区域和形状变化不大, 流速大幅度增加, 但是环流的强度减弱。

本文结论与以上研究结论相类似, 即流场的方向受风场方向的影响较大, 同一区域不同风向会产生不同的流场形态, 空间分布不同的区域在同一风场作用下也有较大的差异, 湖区浅水区域最易受到风场的影响, 一般与风向一致, 湖中心深处往往形成流速较大但与风向相反的补偿流, 最终形成稳定的环流区域。随着吞吐流量增加, 湖区的流速增加, 不会改变环流的位置, 只有局部区域(排干渠口附近)的环流形态发生微小变化, 但由于吞吐量流量的增加, 增强了吞吐流的驱动力, 制约了风驱动力, 导致环流强度减弱。张亚(2014)模拟了4种入库工况下水库的水流特征, 结果表明随着吞吐流量增加至极端量, 流场逐渐由受风生流控制转变为吞吐流控制, 与吞吐流方向一致的区域流速得到加强, 方向不同的区域流速先减小后增加。本研究仅模拟了普通算例吞吐流作用下的水流形态, 并未出现受吞吐流影响较大的区域, 但可以推测, 吞吐流量足够大时, 吞吐流作用会完全取代风生流作用, 极端情况下的流场形态会进行进一步研究。

风驱动力通过改变湖区环流形态, 从而导致污染物在局部区域富集(Ishikawa et al, 2002; 周健等, 2018; Liu et al, 2019)。吞吐流通过水体交换, 进而完成湖泊污染物质的扩散、迁移(姜恒志等, 2013)。蒋鑫艳(2019)等分析了2014—2017年乌梁素海冰封期和非冰封期营养盐的时空分布特征, 得到湖区中部叶绿素a含量对全湖贡献最大的结论, 叶绿素a富集区域与本文A、B区域基本一致; 李兴(2011)等分析了乌梁素海水体总氮、总磷、叶绿素、化学需氧量等的时空分布, 其中排干渠口附近总氮、总磷含量较高, 叶绿素集中分布在湖区中部西大滩和东大滩区, 这些污染物富集区域基本与本文环流区域位置高度吻合。以上研究结果表明本文得到的环流区域的分布是合理的, 从水动力学角度解释了污染物空间分布的规律。乌梁素海风生流及吞吐流对湖区流场影响的研究为污染物富集的成因提供有力的理论支撑, 有关于乌梁素海水质方面的数值模拟会在接下来的研究中进行。

4 结论

本文基于二维水动力模型, 模拟了乌梁素海在定常风作用下湖区的流场形态, 并对不同算例下吞吐流量对湖区流场的差异展开了对比分析, 主要的结论如下:

1) 乌梁素海在无风状态下时, 湖区流速较小, 出现大面积滞流区, 单一的吞吐流无法满足乌梁素海水动力条件, 风驱动力才是主要驱动力;

2) SSW、ENE向2种风场作用下, A区域对应环流形态呈中心对称, B、C区域对应的环流位置、形态结构相似, 但方向相反。浅水湖区最易受到风的影响, 产生与风向相同的水流, 深水区形成与风向相反的补偿流。

3) 吞吐流量的增加会增大湖区整体流速。定常风场作用下, A区域由于受到总排干渠的影响, 环流形态随流量增加发生微小变化, 风生流受到吞吐流影响, 环流强度随流量增加而减弱。

4) 风生流、吞吐流作用下会促使营养盐在局部区域富集, 研究乌梁素海营养盐空间分布时, 需要考虑风生流及吞吐流作用。