台风浪风暴潮作用下三维泥沙数值模拟

朱志夏，熊伟

(1.中设设计集团股份有限公司，江苏 南京 210014; 2.浙江贵仁信息科技股份有限公司，浙江 杭州 310051)

台风引起动力强大的台风浪和风暴潮流很容易造成近岸浅滩大量泥沙的起动，在短时间内造成航道或港池的强淤，使堤防与防波堤等港口航道与海岸及近海工程建筑物遭到严重损坏，同时，风暴增水也会给沿海地区带来安全隐患，因此，对工程海域台风浪、风暴潮以及泥沙输移、港口航道强淤等的数值模拟和预测提出了更高的要求。为了更加准确地模拟台风浪风暴潮作用下的泥沙输移、海床演变、港池航道强淤，有必要进行台风浪风暴潮作用下三维泥沙数值模拟技术的研究。

有关波浪潮流共同作用下三维泥沙数值模拟，已取得了不少成果。朱志夏[1]基于Navier-Stokes方程和质量传输方程，应用多重尺度法、小参数方法及变分原理，导出了波浪作用下的浅水环流方程、不恒定非均匀潮流场中随机波的折绕射联合的数值模式及波流联合作用下的悬沙扩散方程。开发并建立了波流共同作用下二、三维嵌套泥沙数学模型，成功应用于天津新港附近海域疏浚弃土处理的工程。梁丙臣[2]基于国际著名的水动力-生态-悬浮泥沙耦合模式COHERENS-SED和第3代波浪模式SWAN，对COHERENS-SED进行了完善与发展，最终建立了波流联合作用下的三维水动力-悬浮泥沙耦合数学模型，并应用结构化网格二重嵌套的方法，成功应用于黄河三角洲滨海区的潮流和悬沙输运规律的研究。王红[3]采用大、中、小三重结构化网格二、三维嵌套模式，结合第3代波浪模式SWAN和三维水动力-多组分泥沙模型EFDC，将模型应用于潍坊港海域实际工程问题的研究。由于采用结构化网格，工程中关注的防波堤附近水动力的模拟结果有些失真。堵盘军[4]基于第3代波浪模式SWAN和ECOMSED模式的改进和优化，进行了长江口、杭州湾海域三维悬沙数值模拟的初步研究。研究表明波浪在近岸潮滩区域对泥沙启动的作用不可忽视。王平[5]根据考虑流场效应的椭圆型缓坡方程，结合水动力泥沙模型FVCOM(an unstructured grid, finite-volume coastal ocean model)，构建了近岸大范围三维波流耦合模型和波流联合作用下的近岸物质输运模型，分别应用于大连湾海域和旅顺琥珀湾海域，模拟湾内的潮流和物质输运。

综上所述，台风作用下波浪潮流耦合三维泥沙数值模拟技术仍有待进一步研发，例如高分辨率台风风场模拟、波-流耦合计算、细颗粒泥沙的沉降与输移的模拟等等。为此，本文应用中尺度大气模式WRF、第3代波浪模型MIKE21-SW、二维潮流模型MIKE21-HD和三维水动力泥沙模型FVCOM，针对台风作用下泥沙输移、海床演变、港池航道冲淤，通过二次开发，建立了台风作用下波浪潮流耦合三维泥沙数值模拟系统。着重讨论台风作用下三维泥沙输移的模拟，主要解决了波流耦合问题和三维模型中增加了波生近岸流、波浪对泥沙与运动的影响、泥沙的絮凝沉降和高含沙量时的制约沉降机制等影响因素，建立了台风作用下波流耦合作用下三维泥沙数学模型，并将模型应用于连云港港30万t级浅滩深开挖航道工程的研究。

1 三维泥沙数学模型FVCOM改进

1.1 FVCOM泥沙扩散方程

FVCOM采用三维波流耦合作用下的泥沙模型(ROMS_SED)的泥沙计算模式，将泥沙按粒径大小分组，考虑黏性沙、非黏性沙，又称为多组分泥沙模型。其悬沙浓度扩散方程为：

(1)

式中：Ci为第i类泥沙的悬沙浓度；AH为水平涡黏性系数；KH为垂向涡黏性系数；wi为第i类泥沙的沉降速率。

边界条件：在自由表面处不考虑大气沉降的颗粒；在海底边界处，泥沙的沉降和再悬浮过程分别作为水中泥沙的源和汇处理，即：

式中：Di为悬浮泥沙的沉降通量，Ei为海底泥沙侵蚀通量，Ei的表达式为：

式中:Qi为泥沙侵蚀系数；Pb为海底孔隙率；τb为海底剪切率；τci为泥沙起动的临界剪切率。Qi依赖于底部泥沙的物理化学特性，可在实验室中测量得到，根据英国HR实测资料，对于黏性土，Qi=2×10-4～4×10-3kg/m2s，α=1.16。

1.2 泥沙模型的改进

针对连云港港淤泥质海岸开敞海域，模型中需要考虑黏性泥沙的絮凝沉降过程和高含沙量时的制约沉降机制。当泥沙浓度小于制约沉降临界浓度时，絮凝沉降速度计算公式采用彭润泽等基于长江口泥沙絮凝实验的表达式，即：

(2)

当泥沙浓度大于或等于制约沉降临界浓度时，泥沙沉速计算公式为：

(3)

式中：Cgel为胶凝浓度，根据细颗粒泥沙的相关研究，计算中Cgel取为250 kg/m3。

2 三维泥沙模型的建立与验证

2.1 模型建立

数值计算采用二维大、中模型和三维小模型三重嵌套的方法，模拟计算连云港港及附近海域“韦帕”台风作用下波流耦合三维泥沙输移过程[6]。三维模型计算范围(D03)北起岚山头，南至中山河附近，计算域为北纬34°22.1′ ～35°6.6′, 东经116°0.0′～120°8.8′，模型采用非结构化三角形网格，网格最小尺度为20 m，位于连云港主港区和主航道处，网格最大尺度为2 000 m，位于计算域的深水开边界处。三维模型的边界条件由中范围台风作用下二维波浪潮流耦合数学模型提供。

2.2 “韦帕”台风作用下波流耦合三维泥沙模型的验证

为了准确地模拟连云港港及其附近海域“韦帕”台风作用下波流耦合的三维泥沙输移过程，根据天津市气象科学研究所应用WRF中尺度预报模式和美国环境预报中心NCEP历史再分析数据，在天河一号大型计算机上模拟计算了整个“韦帕”台风过程。计算中采用超高分辨率的三重网格嵌套技术，嵌套区域D01、D02、D03的分辨率分别为9、3、1 km，为“韦帕”台风作用下波浪、潮流、泥沙输移的模拟计算提供了准确的驱动风场[6]，计算时间为2007年9月17日0时～9月21日20时，共117 h。

1)三维含沙量验证。

在三维泥沙数值模拟计算中，潮位、潮流、波浪、风场等均由二维嵌套模型提供[6]，三维模型外模计算时间步长取1.0 s，内模计算时间步长选取外模的5倍，垂向平均分为10层。底部摩阻系数范围取0.002 5～0.003 5。根据“韦帕”台风期间实测的波浪、水文资料[7]，进行了“韦帕”台风作用下波流耦合的三维泥沙数学模型的率定与验证。其中1#点位于-3.0 m等深线，2#位于-5.0 m等深线(如图1所示)，每个测站分为表下0.5 m、底上0.5 m共3层进行观测。其中，2#测站的分层含沙量过程的验证结果如图2所示，从验证结果可以看出，计算值与实测值吻合较好。

图1 水文观测点位置Fig.1 Hydrological observation point location

图2 2#观测点含沙量验证Fig.2 Verification of sediment concentration at observation point 2#

受“韦帕”台风影响，2007年9月19日0时“韦帕”台风远在浙江东海海面上，距离连云港海域较远，水体含沙量较小。水体含沙量逐渐增大，到9月20日12点，含沙量增至最大值(如图2、图3所示)。然后逐渐减小，主要特征表现为：最大含沙量时刻落后于最大波高时刻。根据三维泥沙数值模拟结果得到，海州湾内南部的高含沙量水体将随着风暴潮流横向越过主航道后进入徐圩海域，可成为连云港港主航道泥沙回淤的重要来源之一。因此，在“韦帕”台风期间，海州湾内南部的高含沙水体对主航道回淤有比较重要的影响。综上，本文建立的三维泥沙数学模型较好反映了“韦帕”台风期间连云港海域泥沙分布特征及输移过程，可为连云港港30万t级浅滩深开挖航道工程的设计提供技术支撑。

图3 2007年9月20日12时含沙量分布Fig.3 Sediment concentration distribution at 12 o′clock on September 20, 2007

2)航道回淤验证。

根据2007年9月13日(“韦帕”台风前)和2007年9月21日(“韦帕”台风后)连云港港主航道断面水深测量[8]，2次实测的航道水深沿程分布得到(如图4所示)，风前航道水深沿程分布在0～1.5 km和1.5～3.0 km呈阶梯状变化，3.0 km以外水深逐渐增大；由于台风期间航道产生淤积，风后航道水深沿程分布在0～3.0 km内基本一致，向外水深逐渐加大。

图4 “韦帕”台风前后主航道沿程水深Fig.4 Water depth along main channel before and after typhoon Wipha

对比2次水深测量结果可知，在航道里程0～1.3 km内，航道内的淤积强度基本一致，在0.3～0.5 m；航道里程1.3～3.0 km的淤积强度较大，约为前者的2倍。模拟计算得到的“韦帕”台风作用下主航道回淤与实测较为一致，模拟计算的最大实质回淤强度为0.40 m。

3 “韦帕”台风浪作用下三维泥沙数学模型应用

3.1 航道工程方案

连云港港30万t级航道工程由连云港区航道、连云港区外航道、徐圩进港航道组成，其中，连云港区外航道外段宽350 m，底标高-23.0 m，外航道内段宽310 m，底标高-22.5 m；徐圩进港航道宽370 m，底标高-22.0 m，徐圩港区内航道宽210 m，底标高-13.3 m；徐圩二港池航道宽170 m，底标高-11.0 m。徐圩港区防波堤口门分为“平口”、“八字口”，结合不同港内围垦工程方案布置，分为4个方案，本文以方案2、方案4为例进行分析，具体布置如图5所示。

图5 方案2、方案4布置图Fig.5 Layout of plan 2 and 4

3.2 含沙量计算分析

图6～图11给出了工程方案2、方案4连云港海域在“韦帕”台风期间代表时刻的表层和底层含沙量分布情况。2007年9月19日0时，“韦帕”台风远在浙江东海海面上，距离连云港海域甚远，除灌河口含沙量较大外，连云港海域内的含沙量均不大(如图6、图9所示)。随着台风中心逐渐往北移动，台风中心距离连云港越近，台风浪也随之加大，导致近岸水体含沙量逐渐增加，波浪掀沙、波流共同输沙作用明显，9月20日0时表层、底层含沙量均有明显增大(如图7、图10所示)。至9月20日12时，连云港海域含沙量达到最大(如图8、图11所示)，滞后最大台风浪出现时刻约5 h[6]。值得注意的是，由于该时刻的风暴潮流为由北向南的沿岸流，海州湾内的高含沙水体会越过连云港区主航道进入徐圩海域，与实测资料较为相符，即-5.0 m测点的含沙量从最高浓度逐渐减小后，仍然有大幅度增大的趋势。

图11 2007年9月20日12时方案4含沙量分布Fig.11 Plan 4 sediment concentration distribution at 12 o′clock on September 20, 2007

图10 2007年9月20日0时方案4含沙量分布Fig.10 Plan 4 sediment concentration distribution at 0 o′clock on September 20, 2007

图8 2007年9月20日12时方案2含沙量分布Fig.8 Plan 2 sediment concentration distribution at 12 o′clock on September 20, 2007

图7 2007年9月20日0时方案2含沙量分布Fig.7 Plan 2 sediment concentration distribution at 0 o′clock on September 20, 2007

由于工程方案2、方案4港区围填面积的不同，港池内部的水体含沙量也有差异，但与口门外相比，含沙量都较小。一方面，由于各方案中防波堤均已建成，防波堤的防浪效果较好，港池内的波高均不大，口门处波高最大值约为2.2 m，波浪自口门处逐渐递减，港池内最大时刻波高介于1～2 m，和港池外相比小得多，难以引起港池内部泥沙大面积启动和底部高浓度含沙层的产生。另一方面，围填区域面积的不同，造成港池内的环流强度均有不同，其特征表现为港池围填面积越小，环流强度越大。如工程方案2环流强度较大的时刻，港池内的环流区含沙量浓度稍高(图6～8)；工程方案4港区围填全部形成，环流强度最弱，港内含沙量较小(图9～11)。因此，徐圩大环抱防波堤建成后，港池内的水体泥沙主要来源于港外输入和港池内大强度环流造成的浅滩冲刷。

图9 2007年9月19日0时方案4含沙量分布Fig.9 Plan 4 sediment concentration distribution at 0 o′clock on September 19, 2007

图6 2007年9月19日0时方案2含沙量分布Fig.6 Plan 2 sediment concentration distribution at 0 o′clock on September 19, 2007

3.3 回淤计算分析

连云港主航道内外段的沿程回淤强度如图12所示。由图可知，航道最大回淤发生在航道里程3～5 km，相当于旗台北防波堤外1～3 km内(旗台北防波堤距离主航道拐点约2.2 km)，回淤强度约为0.40 m。总体来说，主航道外段回淤小于内段，航道里程5～40 km测点随着距离口门越远，回淤强度迅速减小，然后趋于平缓。由于受防波堤的阻挡作用，旗台港区内部回淤较小，如航道里程0～3 km随着离口门越远回淤越小。

图12 主航道内外段回淤强度Fig.12 Siltation intensity of the main channel

徐圩进港航道方案2、方案4沿程回淤强度如图13所示。可见，徐圩进港航道回淤最大范围出现在航道里程3～7 km，该段航道2个方案回淤强度差别不大；其中，最大回淤位于距离徐圩大环抱口门约5 km处，其风后最大回淤强度约为0.30 m。位于口门内的0、1 km处回淤比口门小，平均回淤约0.15 m。总体来说，徐圩进港航道最大回淤强度比主航道略小，但沿程回淤强度变化较为平缓。

图13 方案2、4徐圩进港航道回淤强度Fig.13 Plan 2、4 siltation intensity of Xuwei entry channel

徐圩一港池航道方案2、方案4沿程回淤强度如图14所示。其中，方案4与方案2相比，港区陆域围填全部完成，港池内水域面积最小，一港池航道回淤强度明显减小，并随着与口门之间距离的增大而减小，尤其是一港池内航道里程0～4 km回淤强度明显小于方案2。

图14 方案2、4徐圩一港池航道回淤强度Fig.14 Plan 2、4 siltation intensity of Xuwei port #1

总之，与徐圩进港航道回淤相比，徐圩港区港池内回淤强度较小，主要是因为在台风期间，防波堤掩护作用较好，港内水体含沙量较小；另外，在台风风速较大期间，港池口门处的潮流流向几乎平行于港池口门，使得直接进入港池内的泥沙较少。虽然徐圩航道靠近南部的灌河口海域，由于“韦帕”台风期间，潮流的主流是由西北向东南沿岸流动，所以，灌河口海域的高含沙水体不易将泥沙输移至徐圩航道工程海域，进一步说明波浪掀沙、波流共同输沙的重要性。

4 结论

1)基于三重网格嵌套的超高分辨率(分辨率分别为9、3、1 km)“韦帕”台风风场，应用第3代波浪模型MIKE21-SW、二维潮流模型MIKE21-HD和三维水动力泥沙模型FVCOM，以二维、三维大、中、小模型三重嵌套的方式，主要解决了波流耦合问题和三维模型中增加了波生近岸流、波浪对泥沙运动、泥沙的絮凝沉降和高含沙量时的制约沉降机制等影响因素，建立了台风作用下波浪潮流耦合三维泥沙数学模型。

2)应用台风作用下波浪潮流耦合三维泥沙数学模型，模拟了2007年“韦帕”台风作用下连云港海域三维泥沙输移过程。2#测站的分层含沙量过程的计算值与实测值吻合较好。“韦帕”台风作用下连云港区主航道的回淤强度与实测较为一致，最大实质回淤强度为0.40 m。较好地反映在“韦帕”台风作用下连云港海域的泥沙运动特征。

3)计算分析了方案2、方案4连云港港海域“韦帕”台风作用下波流耦合三维泥沙运动特征，得到：连云港主航道航道最大回淤发生在航道里程3～5 km，相当于旗台北防波堤外1 km～3km内回淤强度约为0.40 m。受防波堤的阻挡作用，旗台港区内部回淤较小。徐圩进港航道最大回淤值比主航道略小，最大回淤强度约为0.30 m，但其沿程回淤强度变化较为平缓。徐圩港区港池内回淤强度较小，主要是因为台风期间，防波堤掩护作用较好，港内水体含沙量较小。另外，由于“韦帕”台风期间，风暴潮流的主流是由西北向东南沿岸流动，所以，灌河口海域的高含沙水体不易将泥沙输运至徐圩航道工程海域，进一步说明波浪掀沙、波浪潮流共同输沙的重要性。