基于MIKE21模型的椒江口常态下泥沙输运研究

田中仁，张火明，管卫兵，陆萍蓝

(1.中国计量大学 计量测试工程学院，浙江 杭州 310018；2.卫星海洋环境动力学国家重点实验室，浙江 杭州 310012；3.中国计量大学 工程训练中心，杭州 310018)

椒江是浙江地区的第三大河，是比较典型的强潮型山溪性河流。这一区域港口资源丰富，区域经济发达，但其航道受泥沙淤积等问题影响，严重制约区域发展。自1980年开始，椒江区域的研究逐渐得到众多学者的关注。毕敖洪和孙志林[1]对椒江山溪性强潮河口的塑造过程进行了探讨。符宁平和毕敖洪[2]进一步对椒江区域悬沙运动进行了研究。李伯根等[3]基于GIS技术进行数字化冲淤分析，研究了椒江山溪性强潮河口河床冲淤调整过程。夏威夷等[4]探究了该区域泥沙来源和水力活跃度的影响。江晨曦等[5]用Rouse公式估算了椒江口-台州湾海域的悬沙沉降速度。戴玮琦等[6]对椒江河口水沙的特征进行了分析，并对悬沙分布进行了推算。这些研究在一定程度上增加了对椒江入海泥沙运动和河口演变等问题的关注度，并丰富了椒江口-台州湾泥沙输运的研究方法，但未能建立可靠的可视化模型。

本文旨在通过对椒江口-台州湾这一区域常态下泥沙输运状态进行模拟，率定区域数值模型，力求增强模型的区域适用性，为进一步探讨椒江口-台州湾异常天气下泥沙响应提供研究基础。

1 区域简介

椒江区域位于浙江东部温黄平原的北部如图1，椒江上游干流为灵江与下游支流永宁江汇流之后被称椒江，椒江的流域面积约为6 750 km2[7]。椒江在流经海门、牛头颈之后，经由台州湾流向东海。椒江与台州湾河口犹如喇叭形，河口迅速增宽。台州湾水域开阔，与东海有头门、大陈等一众岛屿间隔。椒江该区域主要潮型为不规则半日潮.椒江口台州湾水域的泥沙主要表现为悬沙迁移，其主要源头为河口区域高浓度的悬沙。

本文采用二维水动力MIKE21模型，基于10 m精度的地理数据构建了椒江口-台州湾泥沙输运数值模型，结合地图影像对计算范围内的河道、堤坝、岛屿等进行划分，选取郭聪[8]文献数据2009年4月26日至4月30日数据作为数据来源，以揭示椒江口常态下泥沙输运的时空变化规律，为水利管理部门椒江区域的防洪规划、治理淤泥决策提供科学性依据。

2 模型构建

2.1 研究范围

海洋北部边界为小泽岛一线，东部为东矶列岛-大陈洋外海，南部为黄礁岛-外大陈岛一线，河道边界为海门断面。各水文测验站点如图2。

图2 研究区域及各测点示意图Figure 2 Measurement points in the study area

2.2 网格划分及地形插值

采用1999年台州湾实测水下地形图作为地形资料进行地形插值。其地形图如图3，插值后模型网格如图4。

图3 研究区域地形图Figure 3 Topographic map

图4 网格划分图Figure 4 Mesh plotting

2.3 参数设置

1)时间步长：

模型采用的时间步长为180 s，共1 780步。

2)床面阻力系数：

本文采用较为常见的对模型进行率定的方法，得到该研究区域相对应的Manning数，本文的模型取值为31。

3)水平涡黏系数：

采用Smagorinsky[9]公式计算得到取0.3。

4)初始条件：

初始水位为零，流速为零。

5)计算基面：

采用平均海平面作为计算基准面。

6)边界条件：

①开边界条件，椒江上游边界采用海门潮位站的实测潮位作为控制条件，外部以同一时期的头门、大陈岛和黄琅等测量站的潮位，并适当地由该区域的潮波传播特性做出调整，作为模型外海三条边界的控制条件。

②干湿边界条件，采用常用的干湿判别法[10]，在闭边界上不考虑渗流，假设法向速度与切向流速的法向梯度为零。

水动力模型的参数如表1。

表1 水动力模型参数设置

由于缺乏同时期的实际测量波浪资料，因而本文根据大陈岛往年各年4月份的波浪资料进行统计，得出该区域内的波浪要素，以此作波浪模型的入射条件，来完成计算。采用统计的4月份的波浪要素显示该地区的主要浪向为NE向，出现的频率为33.2%，占到全4月的约1/3。NE浪向的有效波高为1.2 m，其平均周期时5s。次浪向为N向，出现频率为25.5%，其他浪向出现相对较少。

图5 波浪模型有效波高、波向示意图Figure 5 Wave height and direction of wave model

在模拟过程中，采用MIKE21中BW模块，在潮流模型相同计算范围内采用NE方向及有效波高及平均周期作为波浪的入射条件，计算模拟结果作为波浪条件进行输入。

3 水动力过程验证

3.1 潮位验证

本文模拟了2009年4月26日至4月30日大中小潮期间的水动力场，台州湾内的大陈潮位站以及黄琅潮位站的计算值与实测值比较如图6及图7。实测值和潮位模拟值误差一般在25 cm内，可以看到潮位的仿真值与实测值十分吻合，验证了仿真值的可信度。但是，由于采用冷启动，故初始水位稍有偏差。

图6 大陈站潮位验证Figure 6 Tidal level verification at Dachen station

图7 黄琅潮位验证Figure7 Huanglang tidal level validation

3.2 潮流验证

选取了测点t1及测点t3做了流速验证，其仿真结果如图8图9，从中可以看出，流速的仿真值与实测值趋势量值吻合较好，流速的误差控制在10%以内。

图8 t1测点流速验证Figure 8 Flow velocity at t1 point

图9 t3测点流速验证Figure 9 Flow velocity at t3 point

并且还选取了测点t2及测点t3做了流向验证，流向验证结果如图10图11。同样地可以看出流向与实测值趋势较为吻合，除个别时间点外，流向的误差控制在20%以内。可以看出MIKE21模型对于台州湾水动力模型复现性良好，为接下来应用的悬沙输运模型提供了比较准确的水动力模拟条件，同时也印证了该模型对于该研究区域的水动力适用性。

图10 t2测点流向验证Figure 10 Flow direction at t2 point

图11 t3测点流向验证Figure 11 Flow direction at t3point

3.3 水动力模型结果分析

通过模拟研究区域大潮落潮急流矢计算区域动态流场如图12，可以看出此时河口流路和河道的走向较为一致指向口外，椒江口外区域白沙—黄琅一线表现为较为明显的东海域潮流旋转流性质。在动态图中可以看到江山与头门之间受地形影响近似为往复流状态，涨落潮的主流不在一条流路上。由此计算模拟的流场图展示的结果分析可知，MIKE21模型所模拟的流场结果与计算区域的实际状态基本相符。

图12 大潮落急潮位与流矢图Figure 12 Ebb tide flow vector during the tide

4 常态下泥沙输运

4.1 泥沙参数选取

在MIKE21泥沙输运模型中，在液体中泥沙沉速是标识泥沙运动特征的重要物理量，决定着泥沙沉积数量情况，其量值的大小会直接影响泥沙的淤积量大小，因此悬移泥沙沉速的确定是十分重要的。研究水域的泥沙以源于长江杭州湾输出泥沙南下以及海域来沙为主。河流的上游山溪输出泥沙相对影响较小。研究区域内的悬沙主要以淤泥质粉沙，在河流输出淡水时候易于发生絮凝沉降，而且沉速大。在考虑絮凝作用时不适宜直接给定沉速，但能够由沉降系数计算得到实际沉速，此时计算沉速采用5E-4 m/s。见表2。

为描述泥沙在床面上的状态特性(冲刷侵蚀和床层密度)，需要设定床面上的泥沙参数。将冲刷引起的允许最大含沙量设置为50 kg/m3。冲刷设置两个床层，第一床层为软泥层冲刷指数设为8.1，第二床层设定为硬泥层，其冲刷指数设为1.5。

表2 泥沙模型参数

4.2 结果验证

根据泥沙场的参数条件及前一小节水动力流场计算所得出的流速、水深等条件，对台州湾水域悬沙浓度场进行仿真模拟计算，并以江山测点t7点含沙量作为泥沙场验证。图13为验证点泥沙浓度变化。

图13 t7测点泥沙浓度Figure 13 Sediment concentration at t7

比较了验证点含沙量的模拟值与实测值，t7测点含沙量模拟值与实测值误差大多在30%以内，少量时间段的模拟偏差较大，验证点模拟含沙量量级及其变化趋势与实际情况量级大致相同。这也是由研究区域地形数据的选取未能与研究时间段同期地形图及间隔数年地貌可能发生较大变化，以及在床层等参数设置上存在较大差异造成的。但现阶段仿真也足以证明MIKE21数值模型可以较好地仿真出椒江口台州湾区域泥沙分布情况。这也表示基于MIKE21模型的泥沙计算能适应椒江河口泥沙模拟。

5 结 语

本文主要是对模型区域适应性验证，选择了较大的计算区域，首先采用椒江河口台州湾区域实测潮位、流速、流向等参数资料对MIKE21水动力模型进行了率定验证，结果表明数据的计算值与实测值符合较好，能较好的将反映该区域实际测量变化趋势，从而验证了该模型仿真可行性。随后通过对常态下泥沙状态的模拟，经过率定后的MIKE21计算模型对强潮汐河口的数值模拟具有较强的适应性，取得了良好的可视的效果。本文也为接下来对该区域泥沙对台风等异常天气的响应状况的有效模拟及模型应用，提供了理论基础支持。