辽河口红海滩湿地海域潮流及盐度的数值模拟

李晋，乔会婷，徐天平，张洪兴，张明亮，姜恒志

(1.盘锦鸳鸯沟国家级海洋公园管理办公室，辽宁 盘锦 124010；2.大连海洋大学 海洋科技与环境学院，辽宁 大连 116023；3.国家海洋环境监测中心，辽宁 大连 116023)

湿地被称为“地球之肾”，是水生生态系统与陆地生态系统间的过渡区域，能为人类提供大量的食物、原材料及水资源，其在维持生态平衡和生物多样性等方面起着重要的作用。中国河口湿地主要分布在沿海区域，是海洋与河流的过渡地带。由于在河口处存在潮汐、河流径流、波浪等因素相互作用，该水域流态复杂多变。当含有一定盐分的海水随着上涨的潮水上溯到河口内一定距离，便会产生盐楔，即盐水入侵现象[1]。

为掌握河口湿地水域的水动力特性，了解影响湿地植被生长的盐度因子的时空分布情况，国内外学者做了许多相关的研究[2-11]。Feng等[2]针对坡度较缓的湿地建立了一个二维的非线性水流动力学模型(WETFLOW模型)，并用室内模型和野外池塘湿地对其一维和二维形式进行了验证；Abu-Bakar等[3]将TMS (Telemac Modelling System)模型系统应用于英国的塞文河口和布里斯托海峡，研究植被与湿地沉积物间的相互作用对河口环境的影响；Hu等[4]将Delft-3D模型应用于半封闭布列塔尼海湾，探讨了植被对湿地水流的影响；Stark等[5]采用经过校准和验证的Telemac-2D模型，在不同水深和不同植被覆盖度下，对荷兰斯海尔德河口的水动力特征进行模拟；Boumal等[6]收集荷兰西南部的西斯海尔德河口滩涂和沼泽长期的高分辨率水文数据，通过分析发现，无论是在滩涂还是在沼泽植被区，其潮汐幅度与最大流速间存在明显的线性关系，植被区的流速要比滩涂低；Sabokrouhiyeh等[7]利用湿地形状和植被密度的函数关系研究了湿地对废水处理的效率问题；Coletti等[8]采用水-植被-盐度生态水文模型，对澳大利亚西南部一个盐沼湿地进行了模拟，研究发现，水文条件变化与湿地植被保护息息相关。目前，国内在这方面也进行了一系列研究，张明亮等[9]从植被对水流流动特性影响入手，建立了植被作用下的平面二维k-ε双方程数学模型，通过采用拖曳力法和等效阻力系数法处理植被对水流的影响，模拟了刚性植被覆盖的弯曲水槽的水动力特征，得到植被对水位和速度的影响规律；顾峰峰等[10]利用二维数值模型对扎龙湿地进行了模拟，分析了芦苇在湿地蓄、滞洪过程中的作用；Zhang等[11]通过构建植被河道渐变流水动力试验，基于理论分析及试验研究了不同植被参数对植被河道渐变流流场结构、能量耗散、曼宁系数和湍流结构的影响。

辽河是中国七大河流之一，河长约1400 km，流域总面积约22.9万km2，流经内蒙古、河北、吉林、辽宁4个省区，其在辽宁省境内流域面积约6.92万km2。辽河口国家级自然保护区地处辽东湾辽河入海口处，是由携带大量营养物质和泥沙的淡水与海水互相浸淹混合而形成的适宜多种生物繁衍的河口湿地，其湿地面积总量位于世界第二。湿地植物以芦苇、碱蓬植物群落为主，其中，翅碱蓬是河口湿地常见的盐碱地指示植物，具有耐盐、耐旱、适应性广等特征[12]，形成了壮丽的“红海滩”奇景。水分和盐度的消长变化直接影响着翅碱蓬的正常生长。研究发现，翅碱蓬在低盐浓度下生长最好，当盐度含量达到约10时，碱蓬草由绿色变为红色，当盐度大于16时，将会导致翅碱蓬生长受到严重抑制，甚至死亡[13]。对于芦苇来说，主要的影响因素也是盐度和水分，芦苇适合在低盐度的土壤中存活，土壤盐分过高将会抑其生长。近年来，随着湿地污染、海平面上升、海岸侵蚀、湿地资源开发，以及修建拦海大堤、围海造田和建设采油平台等，使得河口水体环境发生改变，并导致翅碱蓬赖以生存的生态环境发生变化，湿地中翅碱蓬大面积死亡，河口湿地部分生态服务功能逐渐丧失，因此，亟须开展辽河口湿地水动力和盐度作用过程的研究，掌握翅碱蓬正常生长所需条件，为辽河口区水资源优化配置及湿地生态功能的恢复提供科学依据。本研究中，基于Mike 21水动力模型，对辽河口红海滩湿地水深、潮流和盐度进行模拟，并采用辽河口湿地区域的水深、速度、盐度实测数据进行了验证。在验证模型的基础上，从数值上分析湿地植被对潮滩区域水动力的影响，并就径流量对湿地水域盐度浓度的影响及盐度的纵向分布规律进行了解析。

1 二维数学模型

Mike 21 是丹麦水利研究所(DHI)开发的一款广泛应用于水动力研究的模型系统，Mike 21 作为一种可以模拟自由表面二维流动问题的模型，可以胜任涉及海洋、河流、河口及内陆湖泊等水体的模拟研究工作[12]。

连续性方程为

(1)

动量方程分别为

(2)

(3)

2 模型的建立与率定

2.1 计算区域岸线及网格划分

模型计算范围为40°18′12″～40°42′38″N、121°01′46″～122°01′53″E的海域，即西至葫芦岛，东至鲅鱼圈连线以北的辽东湾海域，其北端位于湾顶的辽河河口(图1)。模型采用非结构化的三角形网格，网格通过SMS(Surface Water Model System)软件生成，水平方向分辨率灵活可变，对沿岸、河道等地形复杂区域进行针对性的加密。模型计算区域包含9599个网格节点，18 108个三角形单元，模型为显格式，采用的最大时间步长为1800 s，模型在计算过程中根据CFL条件自动调整时间步长。计算区域的网格设置如图2所示，网格可较好地拟合研究区域内复杂的岸线和海底地形。

2.2 边界条件设置

计算域的外海开边界给定潮位，主要考虑了通过TMD(Tide Model Driver)获取的潮位过程线(M2、S2、N2、K2、K1、O1、P1和Q1分潮)。由于研究区域存在潮滩，因此，模型中引入干湿单元判断技术。模型采用冷启动方式，初始水位和速度均设置为0。根据辽宁省水文信息网所提供的历年观测站的径流数据，以及海湾志中辽河口的多年平均径流量，制作流量的时间序列文件作为辽河上游开边界。

图1 研究区域示意图Fig.1 Sketch map of computation domain

图2 计算区域网格水深及鲅鱼圈观测站位分布Fig.2 Mesh, water depth and the measured gauges at computed domain

2.3 模型验证

本研究中，根据锦州港龙栖湾港区起步区填海工程全潮水文观测分析报告中的实测数据对模型进行率定，2个潮位站位H1、H2，4个潮流站位V1、V2、V3和V4的站位位置如表1所示。模型的计算起止时间为2014年5月1日—6月30日。

从图3可知，模型模拟的潮位较实测值来说稍微偏低一些，一方面可能是潮位站位于边界的原因，另一方面可能是TMD提取的开边界潮位过程线不够精确。但无论是振幅还是相位，吻合都是比较好的，误差在允许范围之内。

表1潮位、潮流观测站点坐标

Tab.1Coordinatesoftidallevelandtidalcurrentmonitoringstations

站位station北纬North latitude东经East longitudeH140°47′36.96″121°04′49.98″H240°50′09.06″121°23′41.76″V140°49′13.56″121°08′28.26″V240°48′39.60″121°15′16.68″V340°48′24.00″121°24′20.94″V440°43′50.34″121°23′31.32″

图3 H1和H2站位水位模拟与实测对比图 Fig.3 Validation of the tidal level at stations H1 and H2

辽河口的潮汐类型为非正规半日潮，一天之内有两高两低不等现象，观测海域平均潮差为2.57 m，潮汐强度属中等；两个验潮站平均潮差相差不大，为0.08 m；实测涨潮历时小于落潮历时，涨、落潮平均历时分别为6 h 5 min和6 h 18 min，历时差13 min。可见，模型较准确地模拟了涨落潮过程，以及高(低)潮发生的时刻和潮高。

图4、图5分别为小潮期和大潮期4个潮流站位点流速、流向模拟值与实测值的对比图。由图4、图5可见，大小潮时的流速均有较小的偏差，流向的模拟值与实测值吻合较好。小潮、大潮期间以东部海域V4站位平均流速最大，为0.225 m/s，西部海域V1站位平均流速最小，为0.15 m/s。V1、V2、V3、V4站位平均流速分别为0.147、0.161、0.181、0.225 m/s，总体来讲，观测海域的水流强度为离岸水域(V3、V4站位)强于近岸水域(V1、V2站位)。

图6为辽河口大小潮期间的涨落潮流场矢量图，辽东湾涨潮流的总体流向是东北向。涨潮时，盖州滩东南部和东部潮流的主流路径是从盖州滩东南以东北及北部方向上溯；到达盖州滩尾以东、大辽河河口西侧正对面时，水流方向转为近北向；当继续前进到达盖州滩的东部和东北部时，水流转向北向西北方向上溯，一直到口门形成了从口外至口内的涨潮流。涨潮时，在盖州滩西侧流的方向主要是北向和东北向，由于三角洲地形的影响，外部边坡的流势较弱，产生了大范围的缓流区。大潮期间，盖州滩被两岸潮流所淹没。落潮期间，由于河心滩的发育原因，使落潮流被分为多级分叉水流。其中，辽河口的水流在盖州滩处被分为两股流向海外。偏向盖州浅滩东侧的分叉流为主叉水流，水流流路从河口的南偏东到盖州滩头东侧微转东南向，至盖州滩下部东侧转向南，与大辽河口出口水汇合之后以西南方向泄出外海。靠西侧的分叉流，出流方向是从南向转西南和偏西向，与大凌河口落潮水流交汇。在落潮的初期，虽然盖州浅滩高滩滩面被水淹没，但由于水深较浅，海床对水流的阻力较大，流势弱，从而形成大面积的缓流区。到落潮后期，大范围的盖州滩露出水面，使滩面潮水流入滩两边的冲刷槽，从而加大了槽内的落潮水流动力。

图4 小潮期4个潮流站位(V1～V4)的流速和流向对比图Fig.4 Validation of current speed and flow direction at four measured stations (V1-V4) during neap season

图5 大潮期4个潮流站位(V1～V4)的流速和流向对比图Fig.5 Validation of current speed and flow direction at four measured stations (V1-V4) during spring season

3 红海滩湿地的水动力特征和盐度变化模拟

3.1 数据来源及湿地植被信息的获取

红海滩湿地植被信息的提取是通过遥感技术方法来完成。采用Landsat 8/OLI卫星于2017年6月3日获取的数据作为卫星数据，轨道号为120/32，空间分辨率30 m。通过辐射定标和大气校正对信息进行预处理后，采用遥感常用指数如归一化植被指数、改进的归一化差异水体指数，对典型地物的光谱特征进行分析，构建出用于研究区域地物提取与分类的决策树(图7)，并提取2017年6月辽河口的水生植被信息，提取的植被主要分为芦苇和翅碱蓬两类(图8)。

河口湿地中植被对水体动力的影响由床底摩擦项系数——曼宁系数的变化来表示，曼宁系数取80 m1/3/s，代表整个研究区域的摩擦力。植被区阻力系数取决于水深、植被密度和植被直径，植物作用的曼宁系数nv(与Mike 21中曼宁系数呈倒数关系)计算公式为

(4)

其中：m为植被密度(株/m2)；CD为拖曳力系数；hv为植株的高度(m)；D为植株的直径(m)。

3.2 湿地水动力特征及盐度变化模拟

在验证辽河口水动力模型的基础上，进一步利用曼宁系数nv公式模拟湿地植被作用下盘锦红海滩海域的水流特征及盐度变化规律。芦苇和翅碱蓬两种植被的直径分别为0.6、0.2 cm；植物高度分别为1.5、0.15 m；植物拖曳力系数分别取1.0、0.3；植被密度分别为65、200株/m2。图9为G1和G2两个站位模拟水深与实测结果的对比图，图10为G1站位在有植被情况下的速度实测值与模拟值对比，植物区内速度在涨落潮过程中呈现双峰现象。从图11可以看出，由于湿地植被的存在，消耗了流经植被水流的能量，增加了水流的阻力，导致植被区内的流速明显降低。图12给出了植物区内G1、G2两个站位盐度的变化规律，计算值和实测值较吻合，误差均控制在5%以内。

图6 大小潮期间涨落潮流场矢量图Fig.6 Flow field pattern for ebb and flow tides during spring and neap periods

图7 基于决策树的植被分类图Fig.7 Vegetated classification based on decision tree

3.3 径流量对盐度分布的影响

河流径流量是影响河口处盐度变化的重要因素之一。本研究中通过计算辽河口红海滩湿地区G1、G2两个站位不同径流条件下的盐度大小，分析了径流量变化对湿地盐度分布的影响。选取条件分别为0、30、101、285、450、600 m3/s 的6种径流，具体计算结果见图13。由于G1和G2站位位于潮滩上，当潮位较高时潮水漫滩，潮位较低时处于干滩状态，该时间段图中未有盐度浓度显示。当径流量为0 m3/s时，G1、G2站位盐度值均为32；当径流量为30 m3/s时，G1站位盐度值降为19.96，G2站位盐度值为26.6；当径流量为450 m3/s径流时，G1站位盐度值仅为2，G2站位盐度值为5.7。由此可以看出，不同径流量条件下，盐度的变化规律有明显的不同，盐度的变化随着径流量的增大逐渐减小，径流量越大，平均盐度越小。辽河口红海滩海域上游径流量主要由辽河控制，径流量呈现丰水、枯水、平水状态并呈季节性变化，枯水径流时，上游淡水径流冲淡盐度的作用较弱，盐水入侵距离较大，丰水径流时则正好相反。通过研究径流与盐度变化的关系可知，为保证湿地植被正常生长的盐度需求，防止高盐海水进入翅碱蓬湿地区，在枯水径流时，可采用增加上游下泄流量来控制红海滩湿地的盐水入侵。

图8 植被分布及观察站位分布Fig.8 Distribution of aquatic vegetation and the measured stations

图9 G1和G2站位水深的模拟与实测对比图Fig.9 Comparison of the observed and calculated water depth at locations G1 and G2

图10 G1站位植被作用下的水流速度比较Fig.10 Comparison of current speed at locations G1 with a vegetation

图11 有无植被作用的流场矢量图Fig.11 Flow field pattern in the absence and presence of vegetation

图12 G1和G2站位盐度计算值和实测值的对比Fig.12 Comparison of the simulated and measured salinity values at locations G1 and G2

图13 不同径流条件下G2、G2站位的盐度变化过程线Fig.13 Simulated salinity values at locations G1 and G2 in different runoff conditions

4 结论

本研究中通过利用Mike 21 非结构网格水动力(HD)模块及其附属的盐度输移模块建立了辽河口海域二维水流、盐度数学模型，对西至葫芦岛、东至鲅鱼圈连线以北的辽东湾海域进行了模拟计算，利用外海实测所得潮位和流速资料对建立的二维水流盐度模型进行率定和验证，验证结果较好，表明模型具有一定的模拟能力，且能较为准确地模拟被研究区域的水动力条件。在此基础上，对辽河口红海滩湿地植被区水深、流速、盐度进行模拟验证，并且比较了有无植被分布的流速大小。结果表明，植被的存在能减小湿地水流流速。通过模拟不同径流条件下盐度值的变化发现，随着径流量的不断增加，辽河口海滩湿地海域盐度总体上呈现逐渐减小的趋势，可见上游河流径流量的变化能够引起翅碱蓬植被区盐度的变化，即径流量的变化会直接影响翅碱蓬植被生长。