基于MIKE21模型的大洋河河口水动力与盐度数值模拟研究.docx

《基于MIKE21模型的大洋河河口水动力与盐度数值模拟研究.docx》由会员分享，可在线阅读，更多相关《基于MIKE21模型的大洋河河口水动力与盐度数值模拟研究.docx（15页珍藏版）》请在第一文库网上搜索。

1、基于M1KE21模型的大洋河河口水动力与盐度数值模拟研究贺雪菲（辽宁省河库管理服务中心（辽宁省水文局）IIOOO3）摘要：河口的水动力是一个复杂完整的过程，只有综合考虑径流和潮流的影响方可获得较高精度的模拟成果。以辽宁省大洋河河口为例，构建基于M1KE21的大洋河河口水动力及盐度模型，结果显示具有较高的模拟精度和适应性，可以用于研究区的相关研究。关键词：MIKE21模型；水动力；盐度：大洋河1研究区概况大洋河位于辽宁省东南部，是辽东半岛地区最大的一条独立入海的河流，地理位置为东经122o5124o07,:北纬39049M0049,o其东侧为鸭绿江支流爱河，北以千山山脉分水岭为界与太子河为邻，西

2、侧有大清河、碧流河和英那河，南临黄海。大洋河水系发育良好，水量丰沛，拥有丰富的水资源。大洋河发源于岫岩县偏岭乡境内的一棵树岭南侧，流经岫岩县、东港市，在东沟县的黄土坎入黄海川。流域（包括沙坝河以西区间）的水资源总量31.42亿其中地表水资源量30.73亿人均占有地表水资源量3513亩均占有地表水资源量1489m分别约为全省平均值的5.0倍和3.8倍。平原区地下水可开采量0.33亿m2模型构建2.1 模型及模块选取通过数学模型对河道、湖泊及河口进行数值模拟由来已久，国内外应用最为广泛的水动力及水质数值模拟软件主要包括DH1M1KE系列软件、EFDCQUA12K模型、SMS地表水模拟系统、De1f

3、t3D软件包等。Jeong等基于EFDC模型对GeumRiver进行三维数值模拟分析，模拟计算了不同流速下盐度入侵特征规律及对水体水质的影响程度口。1iu等基于EFDC模型和HSPF模型耦合建立了St.1ouisBay水动力模型，模拟计算水体的水位变化、流速变化、盐度场变化和水温场变化。Zhou等基于EFDC模型建立了Sevem河口拦河坝的水质模型，计算结果显示二维模型和三维模型的模拟盐度分布场基本一致。本研究地点位于辽宁省丹东市，为更好地模拟河口在径流、潮流影响下的水动力情况以及大洋河河口盐度变化情况，采用丹麦水资源及水环境研究所开发的M1KE21模型对河口部分水动力及盐度进行数值模拟。其中

4、，大洋河河口的二维潮流水动力模型采用MIKE21HD(HydrOdynamiCS)模块进行模拟研究工作，大洋河河口的盐度分布模型采用MIKE21AD(Advection-dispersion)模块进行模拟研究工作。2.2 MIKE21HD(Hydrodynamics)模块MIKE21HD(Hydrodynamics)模型采用模型可以基于笛卡尔坐标系或球面坐标系，数值计算方法采用基于非结构网格的有限体积法。其具体计算原理包括(1)基本控制方程、(2)湍流模型方程、(3)底部应力方程、(4)干湿边界的处理。2.2.1 控制方程hhuhv,o+=hStxy22huhuhvu-,.hpn+=fvh-g

5、h-txyxp()xSh2Spsxhx1(sxxsxyy2PoSXPoPoPkSxx+-(hT+-(hT+husSxXyJsShVAhuv6hv-1ha+=-fuh-gh-1-t8XdySyPQSygh2dpyrby1(yxjs2Po力PoPoPOkdydJ+(F)+卷(hyy)+hvs方程中：t为时间；x、y、Z为右手CarteSian坐标系；为水面相对于未扰动水面的高度即通常所说的水位；h为静止水深；u、v、W分别为流速在x、y、Z方向上的分量；Pa为当地大气压：夕为水密度，。为参考水密度；/=2CSin。为COriO1iS力参数(其中C=O.7291（T4sT为地球自转角速率，为地理纬度

6、）；/u和/”为地球自转引起的加速度；、久为辐射应力分量；、葭、为水平粘滞应力项；S为源汇项；（“%匕）为源汇项水流流速。模型求解采用非结构网格中心网格有限体积法求解，其优点为计算速度较快，非结构网格可以拟合复杂地形。1.1.2 湍流模型方程湍流建模采用涡粘理论。亚网格尺度涡粘值由下式给出：A=中?眄瓦其中生是定值，/是特征长度，形变率由下式给出:1u.Suj/IJ2xixiJ,1.1.3 底部应力方程底部应力b=gG遵循二次摩擦定律:其中是阻力系数，%=（%，）是流动速度。滑移速度与底部应力之间的关系为：UTb=JCJ对于二维计算而言，是关于水深的平均速度，阻力系数可以由谢才系数C或者曼宁系

7、数M得到。曼宁数可以由底床糙率长度得到:1.1.4 干湿边界的处理对计算区域内滩地干湿过程，采用水位判别法处理，即当某点水深小于一浅水深Ary（如0.1m）时，令该处流速为零，滩地干出，当该处水深大于62（如0.2m）时，参与计算，潮水上滩。2.3 MIKE21AD（Advection-dispersion）模块盐度的计算遵循对流扩散方程：SSS2S2S五祓）痂+%而其中：S：垂向平均盐度；DxfDv:X和y方向上的扩散系数；2.4 水深及计算网格进行数学模型计算时，依据MIKEC-Map海图数据库及大洋河工程区实测断面地形图确定水下地形，并将水深换算至当地平均海平面。计算网格由三角形单元构成

8、，这样能较好地拟合岸线。依据研究的需要，采用不同的空间分辨率和网格尺度。为节约计算容量，节省时间，对网格进行三层嵌套，工程区外的网格尺度较大，单元格边长约为1.5km5km,对工程区附近的网格进行加密，网格单元格边长约为8001000m。对工程区水域的网格尺度为50300m,以保证精度。所建数学模型网格节点数为10118个，单元总数为18925个。详见图2.4-2研究区域水下地形以及图示所示。图2.4-1研究区域网格方案大洋河卜.地形曲程（BOAbove2416-248-160-8-8-0|-16-8I-24-16-32-24-40-32-48-40-64-72-80-88-72-80-88-

9、Be1ow（a.全局）水下地形高程（m）Above27.525.0-27.522.5-25.020.0-22.517.5-20.015.0-17.512.5-15.010.0-12.57.5-10.05.0-2.5-0.0-2.5-5.0-7.5-Be1ow7.55.02.50.0-2.5-5.0-7.5（b.局部）图2.4-2研窕区域水卜地形以及图示盐度模型范围同二维潮流模型范围一致。2.5 外海边界条件的确定依据MIKE21HD(Hydrodynamics)模型构建的辽宁大洋河河口二维潮流水动力模型西边界为陆边界，东北边界和东南边界为水域开边界，潮流数值计算时需要给定开边界的潮位过程。计算

10、域的潮流场受南黄海旋转潮波系统和东海前进潮波系统共同控制，为了使计算域的潮流场能够反映实际的潮流运动，本模型的开边界根据MIKE21全球潮汐模型提供的调和常数生成外海潮位过程，并结合ChinaTida1的东中国海潮汐数据库提取。边界条件的选取主要考虑水深因素。边界潮位所采用的MIKE21全球潮汐模型使用了ChinaTiaI数据库同化校正，而由于卫星遥感数据在浅水区精度较差，所以全球潮汐模型的精度在浅水区的使用受到一定限制，故本项目中将大范围潮流模型边界定在水深约70m,距离岸线约15Okm处。依据MIKE21AD(AdVeetiOn-dispersion)模型构建的盐度分布模型范围与二维潮流水

11、动力模型范围一致，故盐度分布模型也相应开边界。取上游各级支流输入性的盐度为0,而外海边界的盐度取32%。(即32PSU)。2.6 参数率定盐度模型的糙率系数、涡粘系数、时间步长和干湿边界与水动力的设置相同。区别在于M1KE21AD(AdVeetion-dispersion)模型中的扩散系数采用类比涡粘系数，比例系数取U衰减系数取0。2.6.1 糙率系数计算海域的糙率是个综合影响因素，是数值计算中十分重要的参数，与水深、床面形态、植被条件等因素有关。本模型采用曼宁糙率系数，根据各海域的不同特点，不同区域取值不同。本研究采用附加糙率的概念，与海床的相对起伏变化对应。0.010.01n=0.01+0

12、.02+HH其中，n糙率，M1KE21中M=1mH水深,m在模型率定过程中，根据需要，对manning图进行了调整，见图2.6-1潮流计算所用Manning图。m图2.6-1潮流计算所用Manning图ManningsMm*(13)sAbove9590-9585-9080-8575-8070-7565-7060-6555-6050-5545-5040-4535-4030-3525-30Be1ow25UndefinedVa1ue2.6.2 涡粘系数涡粘系数采用Smagorinsky公式估算，相应Smagorinsky系数取值为0.22m2So2.6.3 时间步长根据模型网格大小、水深条件动态调整

13、模型计算时间步长，使CF1数小于0.8,满足模型稳定的要求，计算时步长在QO1s30s之间。2.6.4 干湿边界对计算区域内滩地干湿过程，采用网格冻结方法处理。当某点水深小于0.005m时，令该网格点为干点，滩地干出，不参与水动力计算；当某点水深大于0.005m但小于0.05m时,令该处流速为零，该网格点仅参与水流连续方程的计算，当该处水深大于0.1m时，该网格点参与计算，潮水上滩。2.6.5 模型计算时间段根据本次规划的要求和附近海域水文实测资料，计算时间从2019年10月15日至11月15日，包括整个大、中、小潮。其中潮流观测时间为2005年10月29日至10月31日。潮位观测时间为200

14、5年10月29日至10月30日。3计算结果与分析3.1 验证资料3.1.1 水动力模型验证资料HO1-I及HO1-2测站的潮流观测时间为2005年10月29日至10月31日，潮位观测时间为2005年10月29日至10月31日。表3.1-1模型验证潮位和潮流测点位置观测时间观测类别测点名称10/29-10/30潮位HO1-I10/29-10/30海流HO1-I10/30-10/31潮位HO1-210/30-10/31海流HO1-2由于盐度分布特性在丰、平即枯水期影响程度不同，因此本模型将模拟各提水泵站、国考断面的盐度分布情况。现有的实测盐度资料为20152018年丰、平、枯及全年1NO34、1No37站点的盐度值、国考断面2017年6月至11月的实测值，并于2023年4月进行实地取样测量，因而本模型的计算时间为2017年1月1日至2018年1月1日，包括整个大、中、小潮。各盐度观测点的坐标如图3.1-2实测站点分布示意图所示。图3.1-2实测站点分布示意图3.2 验证结果3.2.1 水动力模型验证结果(1)潮位验证潮位验证如图321HOI-I潮位验证图所示。潮位验证拟合度较好。HoI-I测站，实测潮位（mHO1T新站：帙拟潮位m图3.2-IH01-1潮位验证图HoI-2测外，实测潮位mHO1-2测站；模拟潮位（m）图3.2-2HO1-2潮位验证图（2）流速及流向过程验证流速及流