临界冲刷切应力对最大浑浊带的影响及其敏感性数值分析

雷　蕾,周志敏,李江文

(1.温州市水乡建设促进中心，浙江　温州　325000；2.温州市水文站，浙江　温州　325000；3.中交水运规划设计院深圳有限公司,广东　深圳　518067)



临界冲刷切应力对最大浑浊带的影响及其敏感性数值分析

雷蕾1,周志敏2,李江文3

(1.温州市水乡建设促进中心，浙江温州325000；2.温州市水文站，浙江温州325000；3.中交水运规划设计院深圳有限公司,广东深圳518067)

受水沙动力和地貌多种因素的影响，河口最大浑浊带内的泥沙处于动态悬浮状态。强切应力场与最大浑浊带相吻合，临界冲刷切应力有规律地先沿程增大而后递减是造成最大浑浊带含沙浓度高的重要原因。为验证临界冲刷切应力对最大浑浊带的位置、分布范围以及悬沙浓度变化有重要的影响，以长江口北槽为例，通过模拟不同临界冲刷切应力下的最大浑浊带悬沙分布，讨论最大浑浊带与临界冲刷切应力的响应关系。研究还进一步表明，滩槽泥沙交换、槽间泥沙交换和底沙再悬浮是北槽最大浑浊带发育与变化的重要原因。

最大浑浊带；临界冲刷切应力；悬沙浓度；数值模拟；长江口北槽

1　问题的提出

由于径流、潮流的相互作用以及盐水入侵等因素的影响，河口区悬沙浓度在纵向分布上往往会出现“上下低，中间泥沙富集”的现象，即存在最大浑浊带[1]。最大浑浊带在河口沉积过程中对泥沙特别是细颗粒泥沙的聚集和沉降起着十分重要的作用，对于潮汐河口而言，最大淤积部位往往取决于最大浑浊带的位置。所以该区段最易出现拦门沙，对入海航道产生明显的影响。

最大浑浊带作为河口地区动力作用最为复杂的区域，许多研究者就其成因、时空变化规律以及对河口淤积的影响等方面开展了大量的工作并取得了一系列有价值的成果。沈焕庭等[2]依据多年的实测资料和数值模拟方法，从径流动力、潮流动力、盐水入侵、悬沙来源和集聚机制等方面对我国一些大型河口的最大浑浊带进行了多角度系统地研究；付德健等[3]通过构建不同形状理想河口的数学模型，研究了河口形状对最大浑浊带的影响；庞重光等[4]探讨了黄河口最大浑浊带的形成机制，指出其形成主要受入海泥沙、泥沙异重流、河口密度环流及湍流的作用；Festa等[5]建立了河口垂向二维数值模型，通过数值计算证实最大浑浊带悬沙浓度位置与水动力条件、泥沙沉降速度、上下游泥沙来源有密切关系；杨云平等[6]基于长江河口1959—2011年实测悬沙浓度数据，研究了长江口最大浑浊带悬沙浓度变化趋势及成因；张文祥等[7]在长江口南槽最大浑浊带地区进行了大、小潮水文泥沙定点连续观测，对长江口南槽九段沙附近最大浑浊带的形成机制进行了初步探讨。研究表明，最大浑浊带悬沙浓度变化具有一定的周期性和规律性。

已有的研究成果表明，河口最大浑浊带内的泥沙处于动态悬浮状态。在该区段存在强切应力场，临界冲刷切应力有规律地先沿程增大而后递减是造成最大浑浊带含沙浓度高的重要原因[8]。临界冲刷切应力通过改变底床的固结程度改变底床泥沙的再悬浮。因此，当改变临界冲刷切应力时，最大浑浊带的悬沙浓度发生改变，进而改变最大浑浊带的时空分布。而目前国内外关于临界冲刷切应力对最大浑浊带影响的研究尚不多见。数值模拟方法作为研究水流泥沙问题的重要手段，已经成熟地运用于最大浑浊带的研究中。但由于泥沙运动过程极其复杂，在构建数学模型时，临界冲刷切应力作为控制泥沙冲刷的关键参数，目前仍需要依靠经验选取[9]，这对数值模拟工作带来不便。若能寻找出临界冲刷切应力与悬沙浓度的响应规律，对于准确模拟最大浑浊带能够提供重要依据。

基于此，本文结合实测资料建立长江口大范围平面二维水流、泥沙数学模型，研究最大浑浊带悬沙浓度变化及发育机制，揭示最大浑浊带与临界冲刷切应力的响应关系，进一步阐明最大浑浊带的形成机理及变化规律。

2　长江口二维水流泥沙数学模型

在长江口深水航道二期工程实施后，位于北槽的深水航道回淤量增大且集中于中段，自北港中下段至南北槽浅滩段形成了纵向延伸距离约25～40 km的最大浑浊带，其含沙浓度表层变化在0.1～0.7 kg/m3，底层变化在l.0～8.0 kg/m3，其最大浑浊带特征明显[10]，历史资料丰富，对开展研究带来了便利。

2.1基本方程

2.1.1水流控制方程

根据Boussinesq 涡黏假定和流体静压假定，建立平面二维水动力数值模型[11]，模型控制方程如下：

连续性方程：

(1)

动量方程：

(2)

(3)

2.1.2泥沙控制方程

悬移质输移扩散方程：

(4)

式中：S为垂线平均含沙量(kg/m3)；Dx、Dy分别为x、y方向的泥沙扩散系数(m2/s)；FS为泥沙源汇函数或床面冲淤函数(kg·m-2s-1)，按式(5)确定：

中国英语学习者的实验在被试所在中学或大学的教师办公室进行，英语母语者分别在各自所在大学的图书馆中进行，一次仅有一个被试在房间中接受测试。被试首先阅读实验要求，然后开始测试。在电脑的自测步速阅读完成后，被试还要做二语水平测试，并填个人语言背景表。二语水平测试题选自Oxford Proficiency Test，共50道语法选择题，用以检测学生的二语语法水平。所有学生均未在之前做过这一测试。语法选择题每题1分，小于30分被界定为低水平;30～35分为低到中等水平;35分以上为中等以上水平。

FS=αω(S-S*)

(5)

式中：S*为水体的挟沙力(kg/m3)，一般采用根据现场资料的经验公式法或半理论方法确定；ω为泥沙沉降速度(m/s)；α为泥沙沉降几率。

2.2模型建立

模型覆盖长江河口区下段至邻近的外海区域：上游西边界至徐六泾，北边界至北纬32.05°，南边界至北纬30.87°，东边界至外海30～35 m等深线处，经度约为122.68°。计算域采用不规则三角形网格剖分，网格在深水航道工程、主槽深泓线处适当加密，最小空间步长为30 m，三角形网格节点36 438个，三角形单元71 406个(见图1)。

图1　计算域网格示意图

数学模型的外海开边界从中国近海潮波模型中提取[12]；上游边界徐六泾，给定2002年2—3月的实测水位变化过程；泥沙边界根据水流饱和挟沙能力公式[13]计算得出。

2.3模型验证

2.3.1水流模型的验证

模型采用2002年2—3月的水文、泥沙实测资料对计算结果进行验证。潮位站和流速测点位置见图2，验证点基本涵括了长江口南支、北支，具有较好的代表性。验证结果表明，数学模型计算结果与实测资料吻合较好，符合精度要求，可以利用此模型开展进一步的研究。部分验证结果见图3～5。

图2　验证测站站点分布图

图3　潮位过程验证图

图4　测站大中小潮流速流向过程验证图

图5　测站悬沙浓度验证图

3　北槽最大浑浊带的特征

3.1北槽最大浑浊带平面特征

图6为模型计算的悬沙分布结果图，直观地反映了北槽最大浑浊带的平面特征，结合已有的研究成果[14]进行分析，北槽枯季悬沙场的平面分布具有以下几个特征：

(1)北槽最大浑浊带位于深水航道工程中下段，其内平均悬沙浓度约为1.0～1.4 kg/ m3；最大浑浊带上游侧落潮悬沙浓度大于涨潮悬沙浓度，下游反之；

(2)在一个潮周期过程中，北槽最大浑浊带的位置随涨落潮的变化而周期性地移动：涨潮时退向口内，落潮时向东南方向延展。在落急、落转涨、涨急、涨转落4个时刻及其附近分别存在悬沙浓度Ⅰ类峰值和Ⅱ类峰值[15]；

(3)在拦门沙段以及进口段，最大浑浊带呈“舌状”分布；在口门外悬沙浓度梯度随水深增大而减小，平均悬沙浓度约为0.20～0.40 kg/ m3。

3.2北槽最大浑浊带与悬沙浓度纵向分布特征

根据北槽2007—2009年悬沙浓度的实测资料，统计分析了北槽12个站点NG3～CS5(见图7)的洪枯季大潮潮周期垂线平均含沙浓度以及潮周期临底平均含沙浓度(见图8)。从图8可以看出，北槽各个站点的悬沙浓度纵向上具有“低—高—低”的分布特点，最大浑浊带位于北槽中段的CS6、CSW、CS3之间。洪季大潮时最大浑浊带内的悬沙浓度高于枯季大潮，洪季大潮的底层悬沙浓度和垂线平均浓度相差较大，底层的悬沙浓度达到了垂向平均浓度的5倍以上(200708的CSW站以及200708的CS6站)。相比较而言，枯季大潮的底层悬沙浓度和垂线平均浓度的差值相对较小，枯季大潮实测最大悬沙浓度介于1.0～1.5 kg/m3，与前面的分析结果相吻合。此外，由于洪季大潮是流域的集中输沙期，输沙量占全年的87%，因此导致最大浑浊带呈现纵向范围向上下游延伸，核心下移的特点[16]。

图7　长江口北槽水文泥沙测站分布示意图

4　临界冲刷切应力与北槽最大浑浊带敏感性分析

4.1临界冲刷切应力场概化

由于受到泥沙粒径、盐度、流速以及水深等因素的影响，临界冲刷切应力是一个非常复杂的参数。部分长江口深水航道泥沙试验的结果[17]表明，长江口临界冲刷切应力在口内约为0.3 N/m2，向外海逐渐增大,其中北槽深水航道治理工程导堤、潜堤以及坝田的临界冲刷切应力大致为3.0 N/m2[18]。笔者依据此就其分布取值情况进行了大量的试验，虽然分布取值后模型计算结果总体较好，但局部结果仍不尽如人意，故本文对临界冲刷切应力场进行了概化：将北槽深水航道治理工程导堤、潜堤以及坝田的临界冲刷切应力取固定值3.0 N/ m2；其他区域取统一值，分4种假定，分别为0.3，0.5，1.0，2.0 N/m2。在其他参数不变的前提下，通过改变临界冲刷切应力，讨论其对北槽最大浑浊带的影响程度。

4.2敏感性分析

为精确分析临界冲刷切应力对长江口北槽最大浑浊带的影响程度，从长兴岛尾部至深水航道末端，沿深水航道的主槽，每隔2 500 m布置一采样点，共计28个点(见图9)。根据枯季大潮的模拟结果，提取4种假定下采样点的悬沙浓度，进行潮周期平均计算，计算结果见图10。

图9　悬沙浓度分析点位置示意图

图10　4种假定下北槽悬沙浓度纵向分布图

由图10可见：假定1(0.3 N/m2)：临界冲刷切应力较小时，床沙较容易悬浮，甚至在水动力较弱的憩流时刻泥沙仍可以较容易地保持悬浮状态，因此悬沙浓度相对较高。北槽最大浑浊带位于距长兴岛尾部40.0～47.5 km段(采样点17～20)，最大浑浊带纵向分布相对集中，潮周期最大悬沙浓度达到1.33 kg/m3；

假定2(0.5 N/m2)：当临界冲刷切应力增大时，悬沙浓度总体下降，最大浑浊带位于距长兴岛尾部37.5～47.5 km段(采样点16～20)，最大浑浊带纵向分布相对集中，但较假定1，位置向上游移动。在最大浑浊带内，潮周期最大悬沙浓度减小至1.17 kg/m3；

假定3(1.0 N/m2)：当临界冲刷切应力持续增大后，悬沙浓度相比前两者下降较多，最大浑浊带的纵向分布范围变广，最大浑浊带位置向上游移动，大致位于距长兴岛尾部32.5～45.0 km段(采样点14～19)，其潮周期悬沙浓度最大已降至0.84 kg/m3；

假定4(2.0 N/m2)：当临界冲刷切应力设定较大值时，床沙相对较难进入悬浮状态，悬沙浓度进一步降低。最大浑浊带的纵向分布范围相比前三者变得更广，最大浑浊带位置向上游移动，大致分布在距长兴岛尾部25.0～45.0 km段(采样点11～19)。最大浑浊带内悬沙潮周期平均浓度最大仅为0.43 kg/m3。

最大浑浊带的变化表现在其含沙浓度、位置和范围的时空变异[19]。从图10可以看出，临界冲刷切应力对于最大浑浊带的模拟结果影响总体较为敏感，尤以最大浑浊带核心部位最为强烈。在不同临界切应力设定下，核心区上游侧悬沙浓度变幅较大，而核心区下游侧至外海的悬沙浓度变化较小，说明在枯季大潮时北槽最大浑浊带的发育和变化过程中，上游侧泥沙输移和底沙再悬浮的贡献大于海相来沙。造成这种现象的原因有2个：一是枯季北槽最大浑浊带具有上游侧落潮悬沙浓度大于涨潮悬沙浓度，下游侧涨潮悬沙浓度大于落潮悬沙浓度的性质。当临界冲刷切应力增大时，床沙再悬浮的动力减弱，北槽整体悬沙浓度减小，上游侧向下输沙和下游侧向上输沙量都减少，削弱了悬沙在高浓度悬沙带的富集，导致核心部位悬沙浓度降幅大于其他部位，最大浑浊带分布范围扩大并且向上游移动，形成了纵向上“低—中—低”的变化趋势；二是南北槽间的九段沙、北槽与北港间的横沙东滩等浅滩发育于汊道间，作为流域来沙堆积的产物和场所，在水动力作用下，其部分沉积物又会发生再悬浮，再搬运和重新分配，成为最大浑浊带次生的补给源。综上说明：枯季大潮时，滩槽泥沙交换、槽间泥沙交换和底沙再悬浮对于长江口北槽最大浑浊带发育与变化的贡献大于海相来沙。

临界冲刷切应力与床沙颗粒组成、盐度、床沙密度以及固结程度等都有直接关系，其应当随位置、水深、床沙组成的变化而改变。虽然，对大范围临界切应力给以固定值的处理方法，会带来一定的模拟误差，但观测资料表明，4种假定中，假定1和假定2较接近实际情况[20]。说明本文的研究仍然较为可靠地反映出临界切应力对北槽最大浑浊带的影响，为今后最大浑浊带的模拟及率定提供较好的经验和思路。

5　结　语

本文建立长江口二维水流、泥沙数学模型，模拟了不同临界冲刷切应力下的最大浑浊带悬沙分布。模拟结果显示：临界冲刷切应力作为判别泥沙冲淤的重要因素，对最大浑浊带的位置、分布范围以及悬沙浓度高低都有重要的影响，尤以最大浑浊带核心部位对于临界冲刷切应力的变化最为敏感。通过对北槽最大浑浊带的形成机理及变化规律的进一步研究还表明：在枯季大潮时期长江口北槽最大浑浊带发育与变化的过程中，滩槽泥沙交换、槽间泥沙交换和底沙再悬浮的贡献大于海相来沙。本项研究对于认识河口最大浑浊带泥沙输运及解决航道拦门沙回淤等问题有积极意义。

[1]沈焕庭,贺松林,茅志昌,等．中国河口最大浑浊带刍议[J].泥沙研究,2001 (1):23-29.

[2]沈焕庭,肖成猷,孙介民.河口最大浑浊带数学模拟研究的进展[J].地球科学进展,1994,9(5):20-25.

[3]沈涣庭.我国河口最大浑浊带研究的新认识[J].地球科学进展,1995,10(2):210-212.

[4]傅德健,朱建荣,沈焕庭.河口形状对最大浑浊带形成的影响[J].华东师范大学学报(自然科学版),2004 (4):72-78.

[5]庞重光,杨作升,张军,等.黄河口最大浑浊带特征及其时空演变[J].黄渤海海洋,2000,18(3):41-46.

[6]Festa J F,Hansen D V.Turbidity maxima in partially mixed estuaries:A two-dimensional numerical model[J].Estuarine and Coastal Marine Science,1978,7(4):347-359.

[7]杨云平,李义天,孙昭华,等.长江口最大浑浊带悬沙浓度变化趋势及成因[J].地理学报,2013 (9):1240-1250.

[8]张文祥,杨世伦,杜景龙,等.长江口南槽最大浑浊带短周期悬沙浓度变化[J].海洋学研究,2008,26(3):25-34.

[9]Lee V,Mineart P,Armstrong L.Calibration of a Sediment Transport Model for San Francisco Bay[J].Estuarine and Coastal Modeling,2004:1007-1026.

[10]I.Brenon ,P.Le Hir.Modelling the Turbidity Maximum in the Seine Estuary (France):Identification of Formation Process[J].Estuarine,Coastal and Shelf Science,1999,49(4):525-544.

[11]Yuanita N,Tingsanchali T.Development of a river delta:a case study of Cimanuk river mouth,Indonesia[J].Hydrological processes,2008,22(18):3785-3801.

[12]沈焕庭,朱慧芳,茅志昌.长江口河口环流及其对悬沙输移的影响[J].海洋与湖沼,1986,17(1):26-35.

[13]DHI.User Guide and Reference Manual of Mike21[R].Denmark:DHI Water & Environment,2005:14-15.

[14]章卫胜.中国近海潮波运动数值模拟[D].南京:河海大学,2005.

[15]顾峰峰,戚定满,牟林,等.长江口悬沙挟沙力公式参数的率定与验证[J].泥沙研究,2010 (5):1-6.

[16]沈焕庭.长江河口物质通量[M].北京:海洋出版社,2001:193-196.

[17]曹慧江.长江口枯季三维流场悬沙数值模拟[D].上海:华东师范大学,2005.

[18]沈焕庭,贺松林,潘定安,等.长江河口最大浑浊带研究[J].地理学报,1992,47(5):472-479.

[19]李九发,时伟荣,沈焕庭.长江河口最大浑浊带的泥沙特性和输移规律[J].地理研究,1994,13(1):51-59.

[20]Thorn M F C,Parsons J G.Erosion of cohesive sediments in estuaries:an engineering guide[C]// Proc.of the Int.Symp.on Dredging Technol.USA:Coastal engineering Proceedings,1980:349-358.

(责任编辑郎忘忧)

Influence of Critical Erosion Shear Stress on Turbidity Maximum Zone and its Sensitivity Numerical Analysis

LEI Lei1，ZHOU Zhi-min2，LI Jiang-wen3

(1.Wenzhou Water Construction Promotion Center,Wenzhou 325000,Zhejiang,China; 2.Wenzhou Hydrological Station,Wenzhou 325000,Zhejiang,China; 3.Shenzhen Branch of CCCC Water Transportation Consultants Co.Ltd.,Shenzhen 518067,Guangdong,China)

Affected by various factors of water and sediment dynamics,and geomorphology,sediment in estuarine turbidity maximum zone keeps in dynamic suspension state.Because the strong shear stress field is in conformity with the turbidity maximum zone,critical erosion shear stress first increases and then descends along the way regularly,which is an important reason for the high suspended sediment concentration in the turbidity maximum zone.To verify the significant impact of critical erosion shear stress on the position of turbidity maximum,distribution range and changes of suspended sediment concentration,north passage of Yangtze River estuary was taken as an example,by simulating the suspended sediment distribution in turbidity maximum zone under different critical erosion shear stresses to discuss the response relationship between turbidity maximum zone and critical erosion shear stress.The study further explained that sediment exchange between channel and shoal,sediment exchange between grooves,and bottom sediment re-suspension are the important reasons for developing and changes of turbidity maximum zone.

turbidity maximum zone;critical erosion shear stress;suspended sediment concentration;numerical simulation;North Passage of Yangtze River estuary

2015-11-12

雷蕾(1988-)，女，助理工程师，硕士，主要从事港口、海岸及近海工程研究工作。

TV131

A

1008-701X(2016)03-0010-06

10.13641/j.cnki.33-1162/tv.2016.03.003