舟山螺头水道悬沙数值模拟研究

2013-08-14冯沈科姚炎明

海洋通报 2013年1期

冯沈科，姚炎明

（浙江大学海洋科学与工程学系，浙江杭州 310058）

螺头水道位于舟山群岛中部，是连接杭州湾和外海的主要潮汐通道之一，其北靠舟山岛南部诸岛，南临穿山半岛，西接金塘水道和册子水道，东连峙头洋，通过舟山群岛东南部诸水道与外海相连，与北侧的长江口和杭州湾存在着复杂的水体和物质交换。螺头水道基本为东西走向，至大榭岛东侧转为西北——东南走向，水深较大，100 m等深线几乎贯穿整个水道，是宁波——舟山港重要的深水航道资源（图1）。

岛礁众多、峡道密布是舟山群岛的一大特点，峡道因其特殊的水动力及沉积特性引起了不少学者的兴趣，对舟山金塘水道、册子水道、虾峙门水道、及马岙水道等的研究已得出舟山峡道水深流急、潮流往复、沉积物具有显著分带性等结论（蒋国俊等，1998），但对于峡道内的悬沙分布及输运特性的研究还相对较少。螺头水道作为舟山岛南侧连接杭州湾与外海唯一的潮汐通道，其内部的悬沙分布及输运特性直接反映了杭州湾与外海之间的泥沙运动和交换特征。已有的有关螺头水道泥沙特性的研究较少，且主要以实测水文泥沙资料分析为主，如有学者分析了螺头水道及其附近海域内部分实测站点的悬沙及底质资料并建立了该水域的潮流模型，从悬沙和底质粒径及水交换的角度探讨了崎头洋附近海域的泥沙特征及悬沙来源，认为杭州湾内的泥沙在水动力的搬运作用下经螺头水道进入崎头洋并向外海扩散（黄惠明等，2009）。然而通过实测资料分析海域内的悬沙特性虽具有一定的合理性，但对资料选取的代表性要求较高，且所得的结论往往存在一定的片面性，因此，本文利用分散在螺头水道内的5个实测站点的水文泥沙资料建立了以螺头水道海域为重点研究区域的二维潮流泥沙数学模型，并在验证良好的模型基础上进一步探讨螺头水道内的悬沙分布及输运特性，加深对于螺头水道内泥沙运动规律的了解，为该海域海岸及航道资源的合理开发利用提供一定的借鉴。

图1 潮流及泥沙调查站位示意图

1 区域潮流模拟

1.1 控制方程

本文采用非恒定流的平面二维有限元水动力数学模型来模拟潮流场运动，此模型被广泛应用到近海潮流数值模拟之中（李孟国等，1999），模型控制方程为：

式中：z为潮位，U，V为x，y方向上的垂线平均流速分量，h为水深，t为时间，c为谢才系数，n为糙率系数，H=h+z；f表示柯氏系数，w为地转角速度，φ为纬度；g为重力加速度，τbx、τby表示x，y方向底床阻力，Ax、Ay则为涡动粘滞系数。

1.2 模型区域及网格化

考虑到模型研究区域螺头水道及杭州湾——舟山海域的地形特点，同时结合后期建立泥沙模型时需要重点考虑的泥沙来源等多方面因素，本次潮流数学模型的计算范围选取较大：长江口上边界取到徐六泾断面；杭州湾内上边界取到澉浦断面；外海开边界中北边界取到与32°16′N纬线重合，东边界与 123°12′E 经线重合，南边界则与 29°16′N 纬线重合。模型计算区域中包含了舟山海域泥沙的主要来源长江口和杭州湾，是后续泥沙数学模型中泥沙输运过程模拟的必要条件。具体位置及范围见图2，模型采用三角形网格进行计算。对模型重点研究区域：螺头水道、舟山海域及杭州湾采用渐变的网格加密处理，整个计算区域共有34 025个三角形单元，18 361个计算节点，计算最小空间步长60 m，计算时间步长取2 s（图3为模型计算网格）。

图2 计算区域示意图

1.3 定解条件

1.3.1 流场初始条件

图3 模型计算网格

陆边界条件：vn=0

1.3.2 开边界条件

模型开边界采用潮位过程控制，根据计算区域附近的相关资料，分析计算8个分潮的调和常数进而预报各开边界的潮位过程来给定开边界条件。

1.4 模型验证

针对螺头水道的地形特点及模型验证的要求，本文收集了2005年冬季螺头水道内5个站位大潮（2005年1月26日15：00至27日18：00）、小潮（2005年1月18日8：00至19日11：00）的潮流悬沙及底沙数据，同时收集了镇海外游山、马目、定海、沈家门、六横、郭巨及公鹅嘴等7个长期和临时潮位站的同步潮位数据作为本次潮流及悬沙模型的验证资料。本次模型计算了2005年1月6日至2月4日一个月720个小时的流场，并采用以上实测资料对模型计算值进行了潮位及流速流向的验证（2#站位于近岸凹湾内，受特殊地形影响流况较复杂，验证时排除，下文含沙量验证时同）。从模拟结果来看，计算结果总体上与实测数据吻合良好，模拟的流场基本能复演研究区域的水动力情况，且可以作为悬沙模型建立的基础。限于篇幅，本文只列出其中两个潮位站及潮流验证站点大潮期的验证结果（图4-5）。

已有的相关研究结果表明，海域中余流的方向和大小一定程度表示着物质输运的方向和程度，因此，本文模拟计算了螺头水道内的余流场（图6），以便与后文悬沙模型中计算的该水道内的悬沙净输移趋势进行对比分析。

图6 螺头水道大、小潮余流分布图

2 悬沙模型建立

2.1 控制方程

本本文采用对流扩散方程模拟海域悬沙的运动，床面泥沙源函数采用切应力法，该模式在河口海岸地区的泥沙运动研究中应用较为广泛（李孟国等，2003），模型的基本控制方程为：

式中：H为流动水深，H=Z-Z0，c为悬沙浓度，Ex、Ey分别为x方向和y方向的悬沙紊动系数，Sc为单位面积上水体与海底的泥沙交换量，t为时间，Me为单位面积上在单位时间内泥沙的冲刷量，τ 为水流对床面的切应力，τ=ρu2*，u*= κVln（H/n），κ为卡门常数，n为有效粗糙高度，u*为底部摩阻流速，ρ为水体密度，τcr为床面泥沙临界起动切应力，τcd为悬沙临界沉降切应力，ωs为悬沙沉速，g为重力加速度。

2.2 定解条件

2.2.1 悬沙场初始条件

2.2.2 边界条件及相关参数确定

陆边界条件，计算时采用法向含沙量梯度为零，即：

开边界条件：

在本次悬沙模型中，作为上游开边界的长江口徐六泾、杭州湾澉浦断面缺乏系统的同步实测资料，通因此需要根据相关理论及历史资料来设计这两个断面的含沙量过程线。

表1为长江口2005年1月及2008年1月的月径流量及月输沙量值（水利部长江水利委员会，2005），由表可知这两个年份长江口冬季的径流量及输沙量值的年际变化不大，均分别为300亿立方米及300万吨左右，因此本文将长江口徐六泾断面2008年1月的实测含沙量过程线（张志林等，2010）经相位调整后作为模型中该断面的泥沙开边界。

钱塘江的径流量及输沙量对杭州湾内的水体含沙量影响较小，因此模型中没有考虑钱塘江径流的作用。本文根据钱塘江河口应用较为广泛成熟的挟沙力公式ψ=kv2/H（林秉南等，1981）以及2003年5月澉浦断面的实测潮流及含沙量值（朱军政等，2010），结合潮流模型计算的该断面的潮流变化来设计含沙量过程线作为模型中澉浦断面的泥沙开边界。

本模型在泥沙计算过程中采用迎流格式，当挟沙水流自计算区域外流入时，采用边界含沙量作为计算条件，若水流自计算区域内流出时，边界节点上的含沙量采用计算结果。

模型中泥沙沉速取0.3×10-4～2×10-4m/s；紊动扩散系数取Ex=Ey=2.5 m2/s；床面泥沙冲刷系数取 Me=1.5×10-4kg/m2·s；卡门常数 κ 取0.4，；有效粗糙高度n取0.000 4 m；泥沙临界起动切应力取τcr=0.2～2 N·m2；泥沙临界沉降切应力取τcd=0.15～1.5 N·m2。

表1 长江口2005年1月及2008年1月月径流量及月输沙量

2.3 模型验证

根据已有潮流流场的计算结果，结合上文所述的二维悬沙对流扩散方程以及经过模型不断调试率定的各项泥沙参数和相关定解条件计算了2005年1月螺头水道及其附近海域的泥沙场，并采用2005年1月螺头水道内的4个潮流泥沙测站同步测量的含沙量数据对模型计算的含沙量结果进行验证（图7-8），由图可知，各测站的含沙量模拟结果总体较好，含沙量的随潮变化趋势与实测值接近，这其中小潮时1#站、3#站及5#站的含沙量验证结果相对较差，由图可见，小潮时，这3个站的含沙量曲线的变化相对大潮来说较为杂乱，这可能是由于小潮时，水体的湍混作用减弱，泥沙梯度增大，进而有利于悬沙的沉降和再悬浮，使得含沙量的起伏变化较大潮明显，而由于二维模型的局限性，本次泥沙模型中难以体现含沙量梯度对悬沙沉降和再悬浮程度的影响，因此小潮时这3个站的含沙量计算值的随潮变化与实测值存在较大的误差。观察验证结果可以发现，观测时期的小潮含沙量大于大潮含沙量，经分析，这主要是由于本次实测站位观测时间的前一个潮周期潮汛明显大于观测时期，大潮汛时较大的潮差可能会在杭州湾内引起大量悬沙悬浮，含沙量增加，而这股高含沙水体通过潮流的反复搬运作用输运到螺头水道需要一定的时间，到达螺头水道时恰逢小潮时期，从而引起了螺头水道内小潮期含沙量异常增大的情况发生。

图7 各测站大潮含沙量验证

图8 各测站小潮含沙量验证

2.4 模拟悬沙场分布特征

通过提取悬沙场的计算结果，绘制了螺头水道海域大、小潮时期涨落急特征时刻的含沙量分布等值线图（图9-10）。由图可见，螺头水道海域由西至东，含沙量逐渐减小，小潮期含沙量大于大潮期；含沙量等值线分布不仅随潮型变化，还和水体运动方向即涨落潮流有关。

大潮涨急时刻，在涨潮流的主流流路上，含沙量等值线向涨潮流方向突出，这主要是因为涨潮流为外海较低浓度含沙水体，主流路径上的水体泥沙含量明显低于同一断面中其他区域的泥沙含量，这也说明了螺头水道内的含沙量大小主要受过境泥沙控制，本地泥沙再悬浮的作用相对较小；大潮落急时刻，与涨急时刻相似，在落潮流的主流流路上，含沙量等值线向落潮流方向突出，亦即主流路径上来自杭州湾方向的相对高含沙落潮流水体运动速度较同一断面上其他水域水体运动速度快，从而增加了这一路径上的水体含沙量，此外，由于落潮流来自相对泥沙浓度较高水域，这使得落急时刻螺头水道内的含沙量等值线向外海移动，即同一位置水域含沙量落急时刻略大于涨急时刻。与大潮时特征时刻的含沙量分布变化不同，小潮时涨落急特征时刻含沙量等值线的变化相对不明显，这主要是由于小潮时水动力条件较弱，水体运动速度慢，涨落潮主流运动速度优势不明显所造成的。

图10 小潮涨落急含沙量分布图

3 悬沙特性分析

3.1 悬沙输运趋势分析

3.1.1 区域悬沙输运趋势

图11为根据悬沙场计算结果绘制的螺头水道海域大、小潮单宽净输沙趋势分布图，由图可知，在螺头水道内存在两股明显相反的泥沙净输移趋势，这两股趋势在水道内南北分明，而这主要是由于螺头水道内的涨落潮主流流路差异所造成的。涨落潮主流流路不同导致在螺头水道内，特别是南北近岸地区，涨落潮流分配极为不均，形成明显的涨潮优势流或落潮优势流，进而使得这些区域内水体出现单向的涨潮流方向输沙或是落潮流方向输沙。在两个净输沙趋势的交汇处，即接近水道中部的深水区域，由于涨落潮流差异相比近岸水域小，因此单宽净输沙量也相对较小，且输沙方向不定，涨落潮方向的净输沙趋势均存在。水道内单宽净输沙趋势的大、小潮差异并不大，尽管小潮时期水道内的水体含沙量高于大潮时期，但是由于小潮时水体的水动力相对较弱，涨落潮优势流没有大潮时明显，因此在水体净输沙量值上较大潮时小。

图11 大、小潮泥沙输运趋势分布图

对比图5所示的螺头水道内的大、小潮余流场分布图，可以发现水道内的水体净输沙趋势和余流分布情况相类似，即各点的净输沙方向基本与该处的余流方向相一致，且净输沙量值大小也和余流强弱程度息息相关，余流较大处净输沙量值大，余流较小处净输沙量值小。因此，在螺头水道内，余流的方向和大小一定程度上表明了物质输运（主要是悬浮泥沙）的趋势和强度。

3.1.2 特征断面悬沙输运趋势

由以上分析可知，螺头水道内涨落潮方向的净输沙趋势均存在，且从数值上看差异并不大，因此需要通过定量的方法判断螺头水道内的悬沙净输运方向，为此在模型中螺头水道近东西口门处分别布置一个断面，根据悬沙场的模拟结果计算了这两个断面的净输沙量，作为分析螺头水道悬沙净输运方向的依据。图12为断面布置示意图，表2为两个断面大、小潮时净输沙量值统计结果（统计结果落潮流方向为正，涨潮流方向为负）。

图12 净输沙量计算断面布置示意图

由表2的计算结果可知，无论大、小潮，断面的净输沙方向均为落潮流方向，即由螺头水道向外海方向输沙。从数值上看，同一潮型条件下，断面的净输沙量值相差不大，这表明在螺头水道内，悬沙的沉降和再悬浮作用并不明显，水道内的海床相对较为稳定，螺头水道主要作为长江口及杭州湾向外海输沙的主要通道而非主要沉积区域；在不同的潮型条件下，大潮的断面净输沙量明显大于小潮，虽然计算的悬沙场中螺头水道海域小潮时的含沙量大于大潮，但由于大潮时的水动力较强，对泥沙净输运量的影响明显大于含沙量的影响，因此大潮时通过螺头水道向外海输运的泥沙量值仍大于小潮。

表2 断面净输沙量计算结果

3.2 含沙量与水动力关系

在海域中，水动力是导致泥沙悬浮输移的最主要动力因素之一，已有的关于海域水体中含沙量与水动力关系的研究成果表明，含沙量的大小通常与潮流流速的大小存在一定的相关性。本文在螺头水道中部由西至东分别提取了4个特征点的大、小潮流速及含沙量值，绘制成同步的过程线，以探求水道内含沙量与水动力之间的关系。由于这4个点分别位于水道内流速较大的中部，水深也较大，涨落潮流分配相对平衡，排除了地形等因素的干扰，因此能基本地反映水道内水动力与含沙量的相关关系。限于篇幅，仅列出大潮时期各特征点流速及含沙量同步过程线图，图中流速涨潮为负，落潮为正。

由图13可见，各特征点含沙量变化周期都约为流速量值变化周期的一倍，这表明在螺头水道内，含沙量与潮流流速存在的一定的相关性。仔细观察还可以发现，各点一个周期内含沙量的最大值和最小值均出现在潮流转流时刻附近1h左右的时段内，即含沙量的变化与涨落潮变化过程存在着较大的关系。由图可知，当潮流由落潮流转为涨潮流时，各点水体含沙量开始从峰值回落直到涨潮流转为落潮流时，含沙量达到最小值，继而随着落潮流过程开始，含沙量逐渐增高，直到下一个转流时刻达到含沙量的峰值。分析原因，主要是由于当螺头水道内水流方向为涨潮流方向，整个海域处于涨潮过程时，涨潮流将外海的低浓度含沙水体带入螺头水道内，使得水道内的含沙量降低，与此过程相反，落潮时期，来自杭州湾方向的相对高浓度含沙水体进入螺头水道内，提高了水道内的泥沙含量。由此可以看出，螺头水道内的悬沙含量受水道上下游水体含沙量变化的影响较大，当地泥沙的沉降和再悬浮程度相对较小，使得含沙量值对潮流流速的变化响应不显著，含沙量值曲线的波动主要反映了螺头水道上下游水体内含沙量对水动力变化的响应，水道内的悬沙输运则以平流作用为主。