大凌河口水动力数值模拟分析

郭常来，崔 健，孙秀波，马玉祥，艾丛芳

1.中国地质调查局 沈阳地质调查中心，辽宁 沈阳 110034，2.大连理工大学/海岸和近海工程国家重点实验室，辽宁 大连 116024

0 引言

大凌河发源于辽宁省凌源市，于锦州凌海市注入渤海辽东湾.大凌河全长447 km，流域面积达2.33×104 km2，是辽西地区最大的河流［1］，其下游及河口三角洲地区的大凌河口湿地，是我国重要的芦苇滨海湿地，具有不可替代的社会经济价值和自然生态价值［2］.近年来，由于大凌河流域内人口增多和工农业飞速发展，以及气候变化等原因，大凌河径流有所减少［3-4］.而在河口地区，由于河流断面扩大，水流速度骤减，因径流减小而产生的问题更为明显［5］.目前，对于河道径流和潮汐的数学模型主要是基于一维、二维水动力数学模型；周跃华等［6］基于GIS与一维水动力模型对陶乐防洪保护区漫溢洪水进行了风险分析；孙玲玲［7］通过建立二维非恒定流数学模型模拟了黄壁庄水库洪水演进过程；锁晓南等［8］基于二维水动力数学模型对黄河四排口河段进行洪水演进模拟分析水动力场变化状况.基于二维水动力数学模型目前广泛应用于河床变形与泥沙运移中，对复杂地形的计算精度高、可靠性强.本研究基于MIKE 21中的水动力模块，考虑大凌河的河道径流和潮汐的影响，对辽东湾北部海域潮汐进行模拟，重点分析大凌河口附近海域的水动力特征，为进一步研究大凌河口海域的污染物输运及泥沙运动提供准确的水动力基础条件（图1）.

图1 大模型计算域地形分布、测站位置及网格划分Fig.1 Terrain distribution，station location and grid division in large-sized model calculation domain

1 研究方法

MIKE 21是一款用于模拟河流、湖泊、河口、海湾、海岸及海洋的水流、波浪、泥沙及环境的工程软件包.本研究采用MIKE 21的水动力模块建立水动力数值模型.该模块基于Boussinesq［9］假定的二维不可压缩的Navier-Stokes方程［10］，其控制方程组如下：

1）连续性方程

2）运动方程

式中：t表示时间；η表示自由水面水位（m）；h表示总水深（m）；ρ表示水体密度；ρ0表示水体相对密度；g为重力加速度；S表示源汇项的流量值；pa表示大气压力；us、vs表示源汇项的流速值；f＝2ωsinφ，其中ω是地球自转角速度，φ是地理纬度；（τsx，τbx）、（τsy，τby）是x、y方向表面风和海底剪切应力的分量；Txx、Txy、Tyx、Tyy是横向应力，包括黏性应力、湍流摩擦、平流摩擦；分别表示x、y方向垂线平均速度，采用（4）式计算.

在该模块中采用基于单元中心的有限体积法进行离散［11］，空间上采用不可重叠的非结构网格划分，能够较好地拟合河口地区复杂的岸线边界，满足对复杂地形模拟计算的需要［12］.

2 研究区域

2.1 模型建立及设置

为了更准确地得到辽东湾北部海域潮流开边界条件，本研究采用大、小嵌套模拟的方法，建立水动力数值模型［13-14］.其中大模型计算域为渤海海域及部分黄海海域，范围大致为117.3—122.4° E、36.5—41.2° N；小模型计算域为辽东湾北部部分海域，范围大致为121.1—122.3° E、40.2—41° N.采用非结构网格对计算域划分，对近岸及地形复杂处进行了网格加密，以更好地拟合复杂的地形条件.在本文中，水深和岸线根据海图确定，同时参考研究地区最新岸线情况进行修正，大、小模型计算域地形分布及网格划分见图1、图2.

图2 小模型计算域地形分布及网格划分Fig.2 Terrain distribution and grid division in small-sized model calculation domain

本模型中设置闭边界处法向速度为零，大模型开边界采用烟台至大连之间的实测连续潮位数据.在大尺度水动力模型验证良好的基础上，导出小模型所需要的潮位边界资料.径流边界主要考虑辽东湾北部沿岸的大凌河、小凌河、辽河和大辽河4条河流的淡水输入，其中小凌河、辽河和大辽河径流采用月平均数据，大凌河径流采用实测日均径流数据［15］.

2.2 模型验证

利用模型对研究海域进行潮位的数值模拟，选取2005年10月大潮期间的模拟值与环渤海8个代表性潮位测站的潮汐表值进行对比，以验证模型的有效性，潮位站点见表1，对比结果如图3所示.可以看出，本模型潮位模拟结果与实测潮汐资料吻合较好，在一个周期内潮位变化基本一致.总体来看，模型能准确地反映渤海海域的水动力变化特性［16-18］.

图3 环渤海代表测点潮位对比图Fig.3 Comparison of tide levels at representative stations around Bohai Sea

表1 环渤海代表潮位站一览表Table 1 Representative tide stations around Bohai Sea

3 讨论与分析

3.1 整体流场特征分析

模型由“冷态”启动，即在数值模型进行模拟计算时，初始时刻计算域范围内水位均设置为0［19-20］，如图4所示.经模拟发现，t＝2 h时水位降低的落潮过程已传播至渤海范围，t＝6 h时落潮过程传播至大凌河河口地区，t＝12 h时落潮后的涨潮过程完全传播至大凌河河口地区.潮位过程每隔6 h涨潮与落潮交替变化，约12 h为一个完整潮位变化周期，每天两次潮起潮落，与实际潮位变化规律一致.

图4 不同时刻流场状态Fig.4 Flow field states at different times

图5给出了模拟期间辽东湾北部海域最低潮位时刻和最高潮位时刻的流场分布结果.可知，辽东湾北部海区最低潮位为－0.82 m，最高潮位为2.13 m.经模拟计算，辽东湾平均潮差为2.4 m.该海区涨潮历时要略短于落潮历时，并且由于河口河道变窄，由外海传入的潮波由于摩擦等作用潮能发生衰减，潮差沿径流向上游逐渐减小.

图5 辽东湾北部海域流场分布Fig.5 Distribution of flow fields in northern Liaodong Bay

3.2 大凌河口潮流分布

通过校验较好的水动力数值模型，模拟了2018年10月辽东湾北部海域的潮流场的分布.模拟结果显示，小潮期海域潮流场分布与大潮期基本一致，小潮期流速较大潮期流速稍小.图6给出了大潮期模拟海域的低潮时、涨急时、高潮时、落急时潮流场分布.

图6 大凌河口海域不同时刻流场分布状态Fig.6 Flow field distribution at different times in the sea area of Dalinghe River estuary

大潮期间低潮时潮流场数值模拟结果表明，大凌河河口附近海域潮流整体流向为南西向，流速普遍低于0.3 m/s，由外海至岸边潮流流速逐渐减小.近岸处表现为沿岸流，流速普遍低于0.1 m/s.此时，在潮流的作用下，大凌河河道径流朝外海加速涌出.

大潮期间涨急时潮流场数值模拟结果表明，大凌河河口附近海域潮流整体由南西向北东流动，速度普遍介于0.3～0.5 m/s之间；小凌河河口至大凌河河口近岸潮流表现为由西至东向的沿岸流，速度普遍低于0.2 m/s.在涨潮的作用下，潮流涌向大凌河河道，大凌河口口门处流速较小.

大潮期间高潮时潮流场数值模拟结果表明，大凌河河口附近海域潮流整体发生倒转，潮流由东向西流动，流速约0.3 m/s.近岸处表现为沿岸流，流速普遍低于0.1 m/s.

大潮期间落急时潮流场数值模拟结果表明，大凌河河口附近海域潮流整体发生转流，整体潮流流向为南西向，流速普遍介于0.1～0.4 m/s之间；在落潮的作用下，大凌河径流加速涌出河道，大凌河口口门处流速约为0.25 m/s.

3.3 潮流动力特性分析

（1）涨潮动力

辽东湾北部海域在涨潮期间潮流流向基本与辽东湾东、西岸线平行，呈现为北东向.由数值模拟的结果可知，在涨潮初期，由于受到大凌河口浅滩地形的影响，涨潮流主要在大凌河口浅滩外侧及西侧近海海域活动.在涨急时刻及涨潮后期阶段，随着涨潮动力的增强，潮流流向主要是北东和北北东向，受河口浅滩地形的影响，潮流流速较弱，形成大范围的缓流区.

（2）落潮动力

在落潮阶段，河口水动力过程主要来自潮水的进退.在落潮初期，河口内浅滩被水淹没，但受地形影响，水深较浅，流速较小，在河口附近形成了大范围的缓流区.到了落潮后期，落潮动力增强，潮间浅滩露出，水流主要经潮沟流出，从而加大了潮沟的落潮水流动力.大凌河的径流量较小，径流流出口门后，河道断面突然展宽，水动力快速消减，难以形成明显的出海主流.

4 结论

基于大凌河口水动力数值模型模拟结果，分析了大凌河口附近海域的水动力变化特征.由于大凌河口紧邻辽河口、大辽河口，该海域属于入海河流冲淡水的交汇混合区，潮流同冲淡水水流相互作用，加之受潮沟密布和潮滩复杂地形影响，形成了该地区独特的流场结构.

大凌河河口附近海域潮流属于辽东湾海流系统的一部分，其主流方向呈北东－南西向，与辽东湾轴向基本一致.涨潮流（北东向）流速略大于落潮流（南西向）流速，最大涨潮流约为0.52 m/s，最大落潮流速约为0.4 m/s，潮流涨落平均潮流强度的分布大体与该海域等深线相适应.