长江口深水航道三维泥沙数值模型相关问题探讨

2020-08-14 11:28顾峰峰

海洋工程 2020年4期

顾峰峰，郭贺

(1. 上海河口海岸科学研究中心河口海岸交通行业重点实验室，上海 201201; 2. 上海海事大学海洋科学与工程学院，上海 201306)

长江口作为长江出海口，受径潮流相互作用，形成了独特的水沙盐输运和地貌特征：中等潮汐强度作用下的外海驱动力，年际变化达6～7倍的上游径流输入，拦门沙区域强劲的盐水入侵、上溯和盐淡水交汇，近底层存在形成机理复杂的细颗粒高浓度泥沙带。长江口历时13年、投资上百亿打造的我国最为繁忙的黄金水道——长江口12.5 m深水航道，这条航道每年花费10～20亿的巨资进行维护，疏浚量曾经达到了1亿m3/a，目前稳定在5×107m3/a左右。在航槽近底处经常能观测到存在约百公斤的高浓度泥沙层，航道回淤的机理十分复杂[1-6]。目前面临着复杂的水动力—物质(水沙)输运机理难以用数学公式描述的巨大挑战，如果有一个三维水动力—物质输运耦合计算的数值模型，能较为全面的揭示其中复杂的耦合作用机理，预测不同工况和假设条件下的主要水文及航道回淤特征，对于开展长江口航道回淤等基础理论研究具有重要意义，一个完善和强壮的数值模型将是研究此类问题的最理想的试验和验证载体。

在过去的二三十年里，河口区域适用的水动力—物质输运耦合数值模型发展取得了巨大的成就。首先，研究的几何空间从平面二维(2D)转向空间三维(3D)，从平面结构网格转向无结构网格；其次，加入了多CPU、GPU计算和机群计算等不断发展的前沿并行计算技术；另外，在求解方法上从Casulli提出的预测步求解等方法[7-10]逐步发展到了可以突破CFL数的欧拉拉格朗日离散方法等等[11-18]。上述的这些数值模型的前沿技术解决或是部分解决了河口区域数值模型计算的“跨尺度”问题[16-18]，即从海洋到近岸的不同尺度带来的计算效率、稳定性、模拟精度、贴合复杂边界等问题。目前比较常用的三维数值模型包括结构网格的POM[19]、TRIM[8]、DELFT3D和ROMS[20]，非结构网格的UnTRIM[21]、ELCIRC，FVCOM[22]、MIKE3、SELFE和SCHISM等；同时国内也有一些类似的开发研究工作[12-15]。上述这些模型都得到了很好的应用，国内几个主要河口的模型应用情况，如表1[23]。

表1 国外模型在国内主要河口的应用情况表Tab. 1 Application of some models in major estuaries in China

虽然数学模型技术已经取得了巨大的成就，并在近年来开展的长江口水沙盐、地貌变化特征及其内在机理的深入研究取得了明显的进展[1-2]，但现有的模型技术在长江口的应用时经常面临应用困难和盐度、泥沙输运模拟精度不高等问题。针对这些问题，很多学者进行了研究，并已取得了不少的研究成果，包括在非结构网格上利用欧拉拉格朗日离散格式、以及并行技术来克服模型计算效率和复杂边界的问题[17-18]。利用有限元方法离散连续方程，减小网格非正交性计算误差，以及对物质输运控制方程利用分步欧拉法和动态时间步长调整的技术手段，解决模型计算的“跨尺度”问题和改善物质守恒性问题。在长江口区域的水沙盐模拟中大量引入盐度斜压和紊动制约的作用机制以提高模拟精度等等[1]。但是，目前这些模型和水沙运动机理的研究在长江口应用中仍然还有不少问题待解决和完善。

1 长江口物质输运模拟计算格式的守恒性和稳定性问题

长江口的物质输运包括盐度、泥沙都是一个长周期的输运过程，通常至少包括了大中小潮的周期性变化过程(约15天)，因而，对于模型计算物质输运的计算效率及守恒性要求较高，另外拦门沙区域存在高浓度泥沙层及垂向高度层化特征又要求计算格式兼具极高的稳定性，这些要求使得模型的计算步长通常较小，从而使得计算量增大和计算效率降低，通常需要引入并行计算的方法来解决长江口物质输运模拟计算中面临的守恒性和稳定性问题。

另一方面，基于水动力的无条件稳定计算格式，研究解决大时间步长条件下的物质输运守恒性和稳定性，也是一个重要的研究方向。从目前的技术条件来看，在无条件稳定离散格式下，如对流项采用欧拉拉格朗日离散格式，可以取较大的计算时间步长，使得水动力计算的效率和稳定性明显提升，但较难满足物质输运的守恒性，通常需要把长时间步长固定或者动态分解成若干个小时间步长后，利用守恒性控制方程进行计算来满足守恒性要求[17-18]。

根据上述，长江口物质输运模拟在计算格式还存在一些问题，但总体上来看，并行化计算的模型开发和无条件稳定的守恒性计算格式，也都已经取得了不少的研究进展，两者各自或者相结合都是研究此类问题的模型发展的重要方向，其研究成果也将大幅提升解决长江口物质输运模拟这类“跨尺度”问题的模拟能力和计算效率。上述问题的示意图如图1。

图1 长江口物质输运模拟计算格式的守恒性和稳定性问题Fig. 1 Conservation and stability of the model scheme for the material transport in the Yangtze estuary

2 河口区域河床切应力的准确计算问题

近底层河床切应力是长江口河床泥沙起悬和落淤的关键参数，其评估方法一般采用对数流速分布曲线来推算近底河床摩阻和切应力(LP方法)[24-29]。然而在河口区域潮汐作用下，采用LP方法的河床切应力评估值的偏差问题十分明显，误差甚至可以达到100%[30-31]，其受非恒定流和垂线层化影响十分明显。长江口区域的河床阻力系数取值通常只有常规取值的1/4，但实际上它并不是一个常数。不准确的河床切应力估算往往导致不准确的河床冲淤变化和近底泥沙通量，同时影响垂线紊流及水动力结构。

实测资料表明，长江口北槽深水航道内强淤积段近底层的局部落急流速甚至可达2 m/s以上，且落潮明显优势，但每年这些区域的航道淤积强度最大可达约10 m/a。当采用常规LP计算方法计算切应力的数学模型模拟计算该区域的河床变化时，往往表现出强烈的冲刷，这充分说明不考虑河口区域河床切应力估算问题的数值模型，其计算结果往往会导致巨大的误差。

图2为常规对数流速分布计算的河床切应力和坐底三脚架基于ADV观测资料推算的河床切应力(TKE法)的比较[1, 24]，两者的平均比值可达4倍多，这一结果说明了两者差异较大。在实际的模型率定和验证时，长江口的三维阻力系数一般取值约为0.000 8～0.001 5，二维计算的糙率取值约为0.011，都明显小于常规取值。图2中TKE法的应力计算公式[24]：

图2 长江口实测河床切应力及常规对数流速分布估算的结果比较(LP和TKE法估算的切应力比较)Fig. 2 Comparison of the bed shear stress calculated by TKE method and by LP method

(1)

分析形成上述差异的原因，主要是由于长江口区域的垂线流速分布在非恒定流、近底密度层化影响及垂线制紊作用等综合作用下偏移了正常对数分布，这一结论符合大多数学者的研究成果[30]。长江口实测流速曲线和基于TKE法的切应力计算得到的常规对数分布曲线[1]的对比如图3，其结果也验证了基于实测紊动切应力推算得到流速曲线相对于常规对数分布曲线已经发生了明显的偏移。

图3 基于ADV实测切应力的长江口实测流速曲线与对数分布曲线对比Fig. 3 Deviation of measured velocity curve and calculated logarithmic distribution curve based on the shear stress measured by ADV in Yangtze estuary

在潮汐非恒定流及紊动制约作用下，针对河床切应力计算方法进行修正，可以修正近底泥沙通量的计算结果，同时也会影响垂向水动力结构和紊动强度的计算，进一步再影响河床切应力和垂线流速偏移，因而这里存在一个循环作用的耦合作用过程。

根据上述可知：水动力—物质输运耦合模型在长江口的应用，需要修正长江口区域的河床切应力估算方法，以解决局部区域动力、近底泥沙通量计算不合理以及多个影响因子之间耦合、循环作用的问题。

一般来说非恒定流条件下的河床切应力与非恒定流的流速过程密切相关，因而非恒定流的切应力修正可以从流速变化过程来入手。为了说明这个问题，根据包含近底层切应力τb的一维非恒定流动量方程，忽略了对流项：

(2)

式中：U为垂向平均流速，η为水面，H为水深，g为重力加速度，t为时间。基于式(2)可得如下式：

(3)

从式(3)可知非恒定流条件下的河床切应力除了和水面坡降有关，还与流速随时间的变化梯度有关，因而河口区域河床切应力计算修正需要考虑径潮流作用下的加减速流的影响，这和很多学者已有的研究结论是一致的。

3 拦门沙区域近底高浓度泥沙形成机理研究及数值模拟问题

拦门沙区域的泥沙问题是长江口三维水沙盐数值模型技术研究中关注的重点问题，长江口近底高浓度泥沙层的存在是拦门沙区域的重要特征(图4)[2]，是导致航道局部超强淤积的重要原因之一，而近底高浓度泥沙层的形成与该区域的盐度斜压力和紊动强度变化密切相关。其主要作用机理目前在理论上可简单描述为：紊动制约导致拦门沙区域的水、沙、盐分层加剧，泥沙、盐度垂线上易于汇聚于近底层，并加强了盐度斜压力作用下水沙上溯能力，即在泥沙垂向汇聚近底层和纵向上溯能力增强两者共同作用下，使得近底层水沙沿航道的纵向向下的输运能力降低，从而纵向上易于汇聚于滞留区域。在文献[1]中，给出了明确的分析结论：拦门沙区域的紊动制约与盐度斜压力共同作用是纵向上形成近底高浓度泥沙的主要原因。

图4 长江口实测的拦门沙段近底高浓度泥沙集中分布形态Fig. 4 The concentrated distribution pattern of high concentration sediment measured in the Yangtze estuary

根据上述，垂向紊动扩散系数是一个描述近底高浓度泥沙形成机理的重要指标参数，其计算模型一般采用通用的计算模型，例如FVCOM和SCHISM等中采用了垂向一维GOTM紊流模型计算包[32]，其包括了常用的双方程模型等。目前不同的紊流模型由于基于不同的假设，得到的计算结果有一定的差异，在长江口的适用性需要详细比较和论证，也需要进行参数的经验率定。不同紊流模型的数值试验的计算结果的差异性，可以参见如图5[33]。

图5 不同紊流模型计算的结果比较Fig. 5 Comparison of calculation results of different turbulence models

为了进一步说明垂向紊动系数的影响，这里选取零方程模型进行实例计算说明，该模型形式简单，在模拟河口物质输运方面有较好的应用；该方程模型中垂向紊动扩散系数和物质输运方程垂向扩散系数仅为Richardson数的函数，其中Richardson数(Ri)定义如下：

Ri=N2/M2

(4)

其中，

(5)

(6)

式中：S和C分别为盐度和含沙浓度；u，v分别为x和y方在垂线位置z处的流速；ρ为浑水密度；ρ0为清水密度。

Richardson数计算式中的M2表示水平流速垂向梯度的作用，N2的物理含义为密度分层的影响，即Richardson数反映了这两种物理作用的相对强弱；一般认为Ri>0.25，即会受到明显的制紊层化作用影响。

动量方程中的垂向紊动扩散系数(Kmv)的计算式为：

(7)

物质输运方程中的垂向紊动扩散系数(Khv)为：

(8)

式中：v0、vb、Kb均为常数，根据Pacanowski和Philander的建议，它们的取值分别为v0=5×10-3、vb=10-4、Kb=10-5m2/s。

根据上述及式(4)～(8)，紊动扩散系数的变化幅度在紊动制约作用下可以减小几百倍。由长江口实测固定点的垂向水文资料推算的紊动扩散系数的计算分布参如图6[1]。根据图6可知，长江口深水航道纵向上在中下段紊动扩散系数出现了明显的大幅度减小。

图6 长江口航道沿程垂线紊动系数(基于实测资料和经验公式推算)和盐度场分布(注：图中实线为盐度等值线，2012年8月17日05:00)Fig. 6 The vertical turbulence coefficient(based on measured data and empirical formula)and salinity along the Yangtze estuary channel (the solid line is salinity isoline, 05:00 17/8/2012)

通常来说，紊动制约作用将大大降低拦门沙区域垂线上水、沙、盐的扩散，这将明显增加垂向上泥沙、盐度的层化现象，并加剧沿航道纵向的盐度斜压上溯作用，两者都将导致航道近底层纵向向外海的输运能力降低，从而促进近底泥沙和盐度更易于汇聚于底部，这又进一步导致垂线密度分层加剧和紊动制约加强，从而形成了一个循环促进的过程和耦合作用的机制，这是形成目前拦门沙河段深水航道近底高浓度泥沙层存在的主要原因之一[1-2, 7]。阐述这一问题的示意如图7所示。

图7 近底高浓度泥沙形成机理示意Fig. 7 Turbulence restriction and formation mechanism of high concentration sediment near the bottom

实际上长江口长期受到风浪的影响(长江口固定水文站——牛皮礁站的年平均有效波高约0.7 m)，风浪条件下的垂向紊动扩散系数将发生明显变化，尤其是航道边滩水深较小使得泥沙垂向分布受到的紊动系数变化的影响加剧，进而影响边滩泥沙输运，主要表现为边滩泥沙紊动加强，泥沙活动性增强。

根据上述，在长江口进行水沙盐的模拟，采用和研发一个能够反映上述近底高浓度泥沙形成机制的模拟技术是十分必要的。

4 长江口近底高浓度泥沙横向入槽模式问题

在长江口人工开挖的航槽内，经常能在横断面上观测到明显的高浓度泥沙大量堆积于主槽的现象(图8)。图8中的近底高浓度泥沙可以在航道内维持相当长的时间，但航道边滩上几乎观测不到。一般认为密度流作用下泥沙归槽输运的影响是重要的成因之一，这一点推论目前尚需要通过更多的数模试验和进一步的现场观测来验证。这里需要补充指出值得关注的有两点：

1) 垂线密度分层导致的紊动制约作用对泥沙横向输运的影响

在中等强度潮汐作用下长江口北槽航道内的垂线紊动强度理论上是较强的。如果没有前述分析得到的结论，即近底层存在很强的制紊作用(紊动系数减小1～2个量级)导致剧烈的水动力及物质密度垂线分层，这样的泥沙汇聚也会瞬时被破坏，无法持续。

2) 近底层高泥沙浓度的制约沉降作用对泥沙横向输运的影响

当近底层泥沙以正常沉速和概率落淤时，理论上分析可知近底层高浓度泥沙形成时间也是较为短暂的，大部分泥沙将在边滩上落淤；然而当汇聚到一定程度时，在形成的高浓度泥沙层内由于制约沉降作用，会使得沉降速度迅速降低(根据室内试验成果可知减小为原来的约1/6[1])，使得更多的近底高浓度泥沙参与到上述由密度差异导致的横向“汇聚”中来，加剧入槽的泥沙量。这一沉降制约的影响同样适用于纵向近底泥沙汇聚的机理分析。

上述泥沙在槽内汇聚的多种影响因子示意图见图8。

图8 长江口深水航道近底高浓度泥沙断面分布实测及横向入槽模式示意Fig. 8 Section distribution and transverse channel model of the highly concentrated sediment near the bottom in the Yangtze estuary deep water channel

长江口北槽航道内的高浓度泥沙在洪季时的实测沿程纵向分布参见图9。从图8和图9及资料分析来看这种槽内高浓度泥沙的存在有长期性和普遍性，这种近底泥沙分布特征是该区域航道局部回淤强度(最大可超过10 m/a)明显大于非航槽区域的主要原因之一。因而图8所示的多种影响因子的作用机制，也是长江口水沙盐三维数值模型及水沙运动机理研究中需要重点考虑和深化研究的重要问题之一。

图9 长江口深水航槽近底高浓度泥沙沿程分布(高低频水深测量)Fig. 9 The high concentration sediment distribution near the bottom along the Yangtze estuary deep water channel measured by high and low frequency sounder

5 结语

基于近年来长江口三维水沙盐数值模型及水沙运动机理研究的发展，系统的梳理了在长江口深水航道数学模型计算中实际所面临的若干个重要问题，主要包括以下四个方面：

1) 长江口三维水沙盐模型离散求解面临“跨尺度”问题，即从海洋到近岸的不同尺度带来的计算效率、稳定性、模拟精度、贴合复杂边界等问题，针对这一问题，需要在现有研究成果的基础上，进一步提升长江口物质输运模拟计算格式的守恒性和稳定性。

2) 近底层河床切应力是长江口河床泥沙起悬和落淤的关键参数，目前河床切应力估算的不准确往往导致河床冲淤变化和近底泥沙通量的不准确，同时影响垂线紊流及水动力结构模拟，因而在长江口应用时，需要修正河床切应力的估算方法，以解决局部区域动力、近底泥沙通量计算不合理以及多个影响因子之间耦合、循环作用的问题。

3) 长江口近底高浓度泥沙层的存在是拦门沙区域重要特征，是导致深水航道局部强淤积的主要原因之一，因而长江口的三维水沙盐数值模型中必须引入近底层高浓度泥沙形成机制，阐明这一机制需要选择和率定一个长江口适用的紊流模型，并需要完善反映紊动制约、盐度斜压力等耦合作用机理和模拟技术。

4) 对长江口深水航道内的高浓度泥沙来说，一般认为密度流作用下泥沙归槽输运的影响是重要的成因之一，这一点推论需要通过理论分析、模拟试验和进一步的现场观测来验证。这里的泥沙归槽运动需要考虑垂向密度分层导致的紊动制约作用和近底层高泥沙浓度的制约沉降作用等耦合作用的影响。