狭长小型河口冲淤特征与演变机制

匡翠萍，赵 钒，顾 杰，汤 俐

(1.同济大学 土木工程学院，上海 200092；2.上海海洋大学 海洋生态与环境学院，上海 201306)

位于陆海交汇处的河口海岸是人类活动最密集的区域.河口海岸的泥沙淤积易引起港口淤塞和生态环境劣化等问题，危害河口地区的健康发展.河口海岸区泥沙输运受径流、潮流和波浪等水动力影响，而泥沙输运引起的海床地形变化及水体中泥沙含量变化又反作用于水流，因此河口海岸区冲淤演变是多动力作用下的综合反应.涨潮(落潮)对应着向陆(向海)输沙，潮不对称输沙是造成潮汐河口海床演变的主要原因[1].近年全球升温、灾害天气频发和高强度人类活动使得河口区水沙运动更趋复杂[2].因此，探明新形势下典型入海河口地貌演变的动力机制受到众多学者关注.相关研究表明，极端天气对海岸的演变和发育有着巨大的影响[3].任美锷等[4]研究了1939年—1981年间风暴潮对江苏省淤积质海岸的影响，指出一次强台风所造成的结果往往超过正常潮汛下整个季节所产生的变化.目前关于开敞型大河口[2,5]的研究较多，而狭长小型河口的研究较少.狭长小型河口和开敞型大河口在地形上有明显差异，导致了水动力差异，从而对泥沙冲淤及演变机制产生不同的影响，尤其是在极端天气(风暴浪)作用下，狭长小型河口的潮流-波浪-泥沙联合作用机制以及该机制对河槽冲淤演变的短时影响值得深入探讨[6].

风河为青岛西海岸新区主要河流之一，风河河口是一个面向开阔海域、口门缩窄、无径流输入且海岸地貌为砂砾质堆积的小型河口.风河河口至上游橡胶坝的河段属感潮河段，受海相来沙的影响，河道淤积严重，限制了河口的综合利用[7].以风河河口为研究对象，采用黄海大模型、青岛西海岸中模型和风河河口小模型三重嵌套 (见图1) 的方法模拟风河河口的潮流场和泥沙场，分析风河河口及其邻近海域的水沙动力特征，探讨风暴浪对无上游径流的狭长小型河口区冲淤演变的短期影响机制.

图1 风河位置和分区、计算区域和大中小模型

1 研究区域概况

研究区域涵盖风河河口及周边海域，西起风河二坝，沿风河入海方向延伸约3 km(见图1)，河口口门宽约300 m.风河两岸经硬化整治，河道顺直，入海口河段基本呈东西走向，入海口上游3.4 km处建有一橡胶坝，橡胶坝以东河段水量不足，淤泥堆积，河口和河道内多浅滩.

研究区域近岸海域受黄海潮波控制，潮汐性质属正规半日潮，最高潮位2.94 m，最低潮位-3.12 m，平均高潮位1.38 m，平均低潮位-1.40 m，最大潮差4.75 m，平均潮差2.78 m(85国家高程基准面).研究区域潮流为正规半日潮，潮流总体表现为顺岸往复流，呈涨潮SW向、落潮NE向的特征.潮流椭圆按逆时针方向旋转，涨潮流速大于落潮流速，大潮期间涨落急流速均大于小潮[8].该海域常浪向为S向，出现频率为47.64%，次常浪向为SSE向，出现频率为12.86%.强浪向为SE向，最大波高7.6 m，次强浪向为SSE向，最大波高6.5 m.此外，在S向、ESE向和SSW向也都曾有4.0 m以上的大浪出现.从各波级的出现频率来看，波高在1.5 m以下的出现频率高达96.6%，大于3.0 m以上的大浪集中出现在SE向—S向.研究区域属砂质海岸，沉积物类型主要有黏土质粉砂、砂-粉砂-黏土混合、粉砂质黏土、中细砂、中粗砂、细砂等八种沉积物，沉积物类型以黏土质粉砂为主[9].

2 模型建立

数值模拟是研究泥沙输运及冲淤演变的一种有效手段[10].采用Mike21软件基于2015年大小潮和设计波况，建立风河河口二维潮流-波浪-泥沙耦合模型.

2.1 模型介绍

Mike21是丹麦水力学研究所(DHI)研发的水环境综合模拟软件，主要模拟河流、湖泊、河口和海岸等水动力环境[11].基于Mike21软件的 FM模块、SW模块和MT模块，建立非结构网格下的双向波流耦合数学模型[12-13].先由FM模块计算得到水位和流速并输入SW模块，后者计算得到波高、波向等波浪参数，同时将辐射应力和底部切应力输入FM模块，得到新的水位和流速代入SW模块进行模拟，两个模块以此方式循环交换数据从而实现双向耦合.波流耦合过程中的辐射应力和底部切应力分别采用二维非线性波理论[14]和综合矢量法[15]进行计算.在这个耦合数学模型中，FM模块和SW模块的时间步长相同，使数据得到同步交换，保证了计算精度.MT模块则通过调用考虑波浪和潮流相互作用机制的波高、水位和流速等参数进行泥沙浓度计算.根据近岸海域泥沙粒径分布特点[9]，研究区域泥沙运动以悬移质为主，海床变形主要由悬移质变化引起.根据泥沙净沉积速率在每个时间步长内对地形网格数据做出更新，更新后的海床厚度输入FM模块和SW模块，从而确保水动力模型的稳定性.

2.2 计算网格、边界条件及参数设定

采用黄海大模型、青岛西海岸中模型和风河河口小模型三重嵌套网格进行计算(见图1).计算网格均采用无结构三角形网格，网格分辨率按研究需求进行控制.黄海大模型北起成山角，南至吕四，向外海延伸约450 km.岸线网格较密，网格水平空间步长为1 240～64 000 m.青岛西海岸中模型北起胶州湾以南，南至风河河口以南30 km，向外海延伸约75 km，计算网格采用无结构三角形网格.对河口区域网格局部加密，网格水平空间步长为144～7 980 m.风河河口小模型西起风河二坝，沿风河入海方向延伸约3 km，对河道内网格加密，网格水平空间步长为27～580 m.综合考虑计算效率和计算稳定性，模型计算时间步长在0.000 1～0.500 0 s范围内自行调节.

黄海大模型共有三个外海开边界(S、E和N)以及岸线闭边界(见图1).外海开边界采用潮位过程控制，潮位过程由Mike21软件包中的全球潮汐预报模型提供，分辨率为0.125°；侧向固边界采用不可滑移条件，即流速为零.大模型给中模型、中模型给小模型提供开边界的潮位和流速流向过程.研究区域风河河口属潮汐河口，落潮时岸滩出露，边界区域采用动边界处理潮间带干湿交换过程，干水深、淹没水深和湿水深分别取0.005 m、0.05 m和0.1 m.由于风河橡胶坝下无径流注入，河流边界采用闭边界.根据相子门站和小麦岛海洋观测站的波浪统计数据[9]，通过频率分析推算外海不同重现期的设计波浪要素.泥沙模型边界由现场实测数据线性差值得到.

模型曼宁数根据计算范围内的粒径资料取平均值70 m1/3·s-1.波浪破碎采用Battjes和Janssen模型，破碎指标选择0.8.床面泥沙干密度取代表性的平均值811 kg·m-3.鉴于研究区域沉积物近岸粒度较粗，由岸向海变细，泥沙中值粒径Φ为0.90～8.60 mm.经过计算和率定调整，设置青岛西海岸中模型的胶州湾以南近岸海域临界冲刷切应力和临界淤积切应力分别为0.45 N·m-2和0.2 N·m-2，其余区域分别为0.72 N·m-2和0.32 N·m-2，冲刷系数取0.000 05 kg·s-1·m-2.

2.3 模型验证

2.3.1潮位验证

黄海大模型潮位验证资料采用2015年潮汐表上青岛测站(见图1)1月7日0时至1月9日0时的潮位过程，测站的位置如图1所示.青岛测站潮位的计算值与潮汐表预报值的比较如图2所示.潮位验证结果表明，潮位的计算值在相位和趋势上与潮汐表理论值误差较小，潮型符合正规半日潮规律，并且落潮历时大于涨潮历时，呈显著潮不对称特征，能较为准确地反映黄海区域的潮位变化过程，可以给中模型提供开边界潮位和流速.

青岛西海岸中模型潮汐验证资料采用文献[16]中S1验证点2009年8月13日11∶00至2009年8月14日11∶00的一次小潮过程，S1验证点的位置如图1所示.S1验证点潮位的计算值与实测值的比较如图2所示，计算值在大小、相位和趋势上与实测值基本吻合.

2.3.2潮流验证

青岛西海岸中模型潮流验证资料采用文献[16]中的S1验证点2009年8月21日9∶00—2009年8月22日10∶00和文献[12]中的S2验证点2009年6月17日17∶00—2009年6月18日13∶00的流速、流向过程，S2验证点的位置如图1所示.流速、流向的计算值与实测值的比较如图3所示，计算值在相位和趋势上与实测值基本吻合.青岛西海岸潮流涨急流速大于落急流速，水深较大的S1验证点较半封闭水域胶州湾内的S2验证点的流速大.

图3 流速、流向、含沙量验证

2.3.3悬沙浓度验证

青岛西海岸中模型泥沙验证资料采用S1验证点2009年8月13日11∶00至2009年8月14日11∶00和文献[12]中S3验证点2000年2月13日12∶00至2000年2月13日13∶00的悬沙浓度(质量浓度)实测资料，S3验证点的位置如图1所示.悬沙浓度的计算值与实测值的比较如图3所示.泥沙本身具有复杂的物理、化学和水力特性，在海水中更有别于在淡水中的属性，目前对它的物理机制还没有完全了解，对泥沙的研究仍处于半经验半理论阶段，因此用数值方法来模拟悬沙浓度过程相对较困难.从悬沙浓度的数值上看，计算值与实测值的量级一致，量级在1×10-2kg·m-3.

为了确保数学模型计算结果能够准确用于研究，模型必须满足稳定性、收敛性和精度要求.采用误差指数法进行模型计算效率评价.百分比偏差(PBIAS)法是误差指数法的一种，表示计算模拟结果与实测值之间的百分比偏差率.百分比偏差法的计算式如下所示：

(1)

式中：M为实测值;S为模型计算结果.根据计算的PBIAS值可将模型计算效率分为优秀(<10)、非常好(10～20)、好(20～40)和差(>40)四个等级.此法用于评价泥沙模型时评定标准参考值为55，当泥沙模型计算模拟结果与实测值之间的百分比偏差率小于55时，认为模型达到精度要求[17],结果如表1所示.

表1 PBIAS法评价

由表1可以看出，潮流模型的PBIAS评价均为好、优秀、非常好，泥沙模型的PBIAS评价均为好和非常好，说明建立的黄海大模型和青岛西海岸中模型是合理的，能够模拟该区域的水沙动力特性.

3 讨论与分析

为探究潮流-波浪-泥沙联合作用机制以及该机制对河口区演变的影响，将风河河口划分为A、B两个区域(见图1)，其中A区为风河入海区，B区为河口通道，重点分析这两个区域的水沙动力和冲淤演变特征.

3.1 潮流特征

图4为风河河口大潮典型时刻流场.受风河地形岸线的影响，河口段潮流沿河道走向呈东西流向，涨潮流向西，落潮流向东，主槽涨急流速大于落急流速.大潮涨急时刻，A区水流受外海潮流动力影响和主槽束水作用，以较大流速流向B区，A区北岸入海口岸线光滑，凸岸挑流效应弱于南岸，因此主流在A区南岸流速骤增至1.11 m·s-1，B区水流较弱，潮流流速约为0～0.32 m·s-1.大潮落急时刻，B区河段较弱的水流向东流入A区，在入海口南岸由于凸岸挑流作用流速骤增至0.80 m·s-1.小潮期间涨落急流速均小于大潮，并且小潮期间B区几乎无潮流流入.在风河河口内和外海侧灵山湾附近分别取点FH和LSW(位置见图1)进行分析.由图5中FH处和LSW处的潮位过程对比可以发现，与LSW处潮位过程相比，FH处潮波振幅减小，而落潮历时增大，即潮波从外海向风河河槽传播过程中，涨落潮历时不对称性加剧，符合在半封闭水域中，潮波经缩窄口门(峡道)后“振幅减小，相位滞后”的特征，即峡道效应.潮波从风河A区向上游B区传播过程中，由于潮汐作用减小并且风河无径流注入，潮流能力逐渐消耗，流速减小至零.B区水深较浅，岸滩时常出水，越往上游，潮汐引起的断面面积变化越大，涨落潮历时不等越严重.

a 涨急

b 落急

图5 风河河口内外潮位计算过程

3.2 河床冲淤演变特征

3.2.1纯潮流作用

图6为大潮期间涨落急时刻的悬沙浓度分布.大潮期间A区南岸受凸岸挑流作用，较强的潮流冲刷河槽，潮流挟沙能力显著，涨潮期间悬沙浓度最高可达0.150 kg·m-3,落潮期间最高为0.064 kg·m-3.B区涨落潮流速小于A区，潮流挟沙能力显著降低，涨潮期间悬沙浓度最高为0.005 kg·m-3,落潮期间悬沙浓度几乎为零.小潮期间河道内的悬沙浓度总体都很小，A区涨潮期间悬沙浓度最高为0.069 kg·m-3，落潮期间悬沙浓度几乎为零，B区涨落潮期间悬沙浓度都几乎为零.

图7为大潮48 h河床冲淤分布图.大潮期间潮流沿岸输沙使得风河近岸海域呈淤积态势，由于挟沙能力在河道各段存在差异，河道内冲淤结果不一，基本呈“上淤下冲”的演变态势.A区冲刷量约为281.58 m3，冲刷深度自南向北减小，南岸冲刷深度为0.015 m.B区淤积量约为513.49 m3，平均淤积厚度量级为10-4m，自西向东淤积高度增加，在靠近A区处最厚可达0.001 m.小潮期间总体冲淤分布特征与大潮类似，但是冲淤量级基本在10-4m及以下.

风河属典型狭长小型潮汐河口，在纯潮流作用下，河口内输沙受到涨落潮流速和涨落潮历时共同控制，采用Brown等[18]提出的断面潮不对称比的方式判断风河河口断面输沙方向.由于B区仅在大潮期间、潮位高于潮滩时才有潮流流入，因此基于大潮期间的潮流泥沙数据，选取四个典型断面C1～C4(见图1)，即B区西侧断面C1、B区弯段断面C2、B区东侧断面C3和A区断面C4，分别计算各断面的涨落潮最大流速比Up=(Vp-flood/Vp-ebb)和涨落潮历时比Td=(Tflood/Tebb),其中Vp-flood和Vp-ebb分别代表涨、落潮期间最大流速，Tflood和Tebb分别代表涨、落潮期间流速超过泥沙起动流速的时长.在模型计算中，临界冲刷切应力控制着泥沙起动与悬移输运，因此Tflood和Tebb为涨、落潮期间底床切应力大于临界冲刷应力的时长.各典型断面Up和Td的关系如图8所示.Brown等[18]总结了输沙不对称的标准：Up<1.0(在象限1中Up<1.2)时输沙方向为落潮方向，即象限3、象限4和象限1黑色实线以下区域均代表落潮输沙占优.图8中，在Td=1线上的点(断面C1～C3)表示涨潮和落潮过程中底部切应力一直小于临界起动切应力.风河河口典型断面均落在象限1和象限2，表示风河河口为涨潮向输沙.在A区断面C4，落潮期间边滩出露，降低了断面落潮挟沙能力，导致在潮周期内输沙主要存在于涨潮过程.泥沙在涨潮过程中涌入上游B区，B区断面C1～C3在潮流过程中受底摩擦和浅水作用表现为涨潮占优(Up>2.0)，并且底部切应力在B区始终小于临界冲刷应力，导致该段河道在涨落潮期间均存在淤积.因此，风河河口区淤积主要受潮流控制，潮不对称输沙是造成河口地貌演变的主要原因.该规律也可体现在其他无上游径流的狭长小型河口，如漳卫新河河口[19].

a 涨急

图6 风河河口大潮期间涨落急时刻悬沙浓度场

Fig.6Distribution of suspended sediment concentration at maximum flood and ebb during spring tide for Feng River Estuary

Fig.7 48 h bed evolution during spring tide (positive values represent deposition and negative values represent erosion)

图8 不同断面最大流速比与涨落潮历时比

Fig.8 Peak flow velocity ratio versus duration ratio at different sections

3.2.2波流耦合作用

极端天气下的波流非线性关系显著增强是导致深槽短时间剧烈变化的主要原因.为分析不同波况对风河河床冲淤演变的影响，以2015年2月大潮潮流过程为基础，通过频率分析推算外海不同重现期的设计波浪要素(波向、波高H13%和波周期T)，设计了如表2所示的四个波况，重点讨论50年一遇设计波浪要素耦合潮流作用对河床演变的影响.由表2可知，受到岸线走向的影响，NE向波浪的风区较短，波高相对较小，而与岸线走向夹角最大的SE向波浪波高最大.

表2 不同波浪要素对48 h河床冲淤变化的影响

由模拟结果可得，波流耦合作用下河口区泥沙冲淤分布与纯潮流作用下大致相同，但在近岸海域两者差异较大.进一步提取不同波况下河床A区和B区的泥沙冲淤变化量和变化率(与纯潮流相比)，量化不同波浪要素对河床冲淤的影响.由表2可知，波况1作用下，A区和B区泥沙冲淤量变化率均小于15%.波况2～4耦合潮流作用下A区冲淤量变化率均小于18%，但波浪作用对B区河床冲淤量影响较大，变化幅度均大于50%，其中波况2和波况3作用下冲淤量变化率超过200%.

在波况1作用下，与纯潮流作用相比，河口区冲淤量变化较小，A区的侵蚀量减轻1.3%，B区淤积量增加14.0%；风河河口近岸海域波浪掀沙效应显著，距岸线垂直距离1 km范围内有大面积冲刷，冲刷深度为0.008 m，口门北侧出现局部淤积，淤积厚度约为0.003 m.

在波况2作用下，与纯潮流作用相比，A区的侵蚀量增加14.1%，B区冲淤量增加约208.1%；波浪在近岸海域掀沙效果显著，海床高度变化在-0.033～0.095 m，口门近海侧出现较强淤积，淤积厚度约为0.061 m.

在波况3作用下，与纯潮流作用相比，A区的侵蚀量增加17.9%，B区冲淤量增加约274.2%；近岸海域大部分区域呈冲刷状态，海床高度变化在-0.100～0.031 m，口门北侧出现局部淤积，淤积厚度约0.003 m，口门南侧出现局部淤积，淤积厚度约0.015 m.

在波况4作用下，与纯潮流作用相比，A区的侵蚀量减轻14.6%，B区冲淤量增加56.5%；近岸海域大部分区域呈冲刷状态，海床高度变化在-0.044～0.012 m，口门北侧出现局部淤积，淤积厚度约0.002 m，口门南侧出现局部淤积，淤积厚度约0.009 m.

波浪在近岸区域引起了大范围的冲刷，48 h海床高度变化的量级达到10-2m.由于波浪作用在口门处几乎削减为零，口门外被掀起的泥沙大部分落淤在口门附近，小部分泥沙在潮流挟带作用下进入口门后逐渐落淤，因此风河河口区在50年一遇波浪作用下并且波浪传播方向与岸线接近垂直时冲淤显著.在三个方向50年一遇波浪中，SE向和E向波浪对风河河床冲淤影响最大.B区淤积厚度最高比纯潮流作用下增加一个量级.

为探究波流耦合作用下狭长小型河口冲淤原因，进一步分析风河河口区48 h冲淤过程中潮流及悬沙变化过程.图9给出了不同波浪要素下A区和B区两个分析点(位置见图1)流速及悬沙浓度变化.从图9可以看出，波浪作用下，涨潮流流速显著提高，落潮流流速变化较小，潮不对称加剧.波况1作用下，A区涨急时刻流速较纯潮流增加约15%，悬沙浓度变化较小，B区流速和悬沙浓度略有减小.波况2作用下，掀沙能力显著强于波况1，A区涨急时刻流速较纯潮流增加约40%，悬沙浓度增加约200%，B区涨急时刻流速较纯潮流增加约50%，悬沙浓度增加约50%.波况3作用下，波浪掀沙效果显著，A区涨急时刻流速较纯潮流增加约25%，悬沙浓度增加约360%，B区涨急时刻流速较纯潮流增加约50%，悬沙浓度变化较小.几乎与岸线平行的波况4作用下，A区和B区潮流流速未见显著变化，但悬沙浓度均显著提高，A区悬沙浓度提高约100%，B区悬沙浓度提高约25%.

a T1(A区)

b T2(B区)

图9 不同波浪要素下A区和B区分析点流速及悬沙过程

Fig.9 Tidal current velocity magnitude and sediment concentration at analysis points in A and B areas under different wave conditions

综合表2中不同波浪要素下A区和B区的冲淤结果发现，狭长小型河口由于口门走向的特殊性，NE向波浪对河口内潮流影响较小，使得口门外掀沙难以输运至河口区内，SE向和E向波浪均可通过促进风河河口区水动力增强来提高潮流输沙能力，故在表2中波况4B区冲淤变化量显著低于波况2和波况3.与纯潮流作用下相比，强浪(波况1)虽然使得A区水动力略有增强，但是就整个河道而言，悬沙浓度变化较小，故在表2中波况1作用下，A区和B区泥沙冲淤变化率均小于15%，即2年一遇强浪耦合大潮作用下，48 h河口区泥沙冲淤主要受潮流动力作用.在风暴浪(波况2～4)作用下，随着近岸底部切应力的增加，口门外近岸海域波浪掀沙作用显著增强，河口口门处悬沙浓度提高了约100%～360%，48 h淤积变化量最高可达274.2%.风河河口对短期极端天气(风暴浪)有敏感的冲淤响应，这也是风河河口区发育有水下浅滩、砂坝和席状砂等典型浪控沉积地貌的原因.

综上所述，狭长小型河口冲淤在正常天气下受潮流控制，而极端天气下受风暴浪控制.由于狭长小型河口口门走向的特殊性，当波浪抵达岸线的入射角与岸线接近垂直时(SE向和E向)对风河河床冲淤的短期影响最大.该规律也可体现在其他狭长型海湾，如唐岛湾[27].风河河口区对大潮期间SE向和E向波浪有敏感的冲淤响应，未来的河口综合整治过程中，应考虑清淤、边滩开挖、导堤建设等治理工程对潮不对称调整的复杂影响，极端天气如风暴浪期间的工程效应也值得考虑.

4 结论

(1) 在河口地形岸线的影响下，河口区潮流特征差异较大.入海口处河道由外海骤然缩窄，流速较外海增大，最大流速可达1.2 m·s-1.河口区内属浅滩地形，仅在潮位高于潮滩时才有潮流流入，使得河道上游流速迅速减小至零.

(2) 纯潮流作用下，水流挟沙能力在河道各段存在差异，总体呈现“上淤下冲”的演变态势，入海口处水流流速较大，为主要冲刷区域，冲刷深度自南向北减小；河口通道水流较弱则为主要淤积区域，淤积厚度自西向东淤积高度增加.

(3) 风河河口冲淤演变在正常天气下受潮流控制，潮不对称输沙是造成河口地貌演变的主要原因，而在极端天气下受风暴浪控制.由于狭长小型河口口门走向的特殊性，抵达岸线的入射角与岸线接近垂直(SE向和E向)的波浪对风河河床冲淤的短期影响最大.

(4) 常浪与潮流耦合作用下，河道内冲淤分布与纯潮流作用下相似.50年一遇波浪耦合潮流作用下，冲淤演变显著变化；近岸海域主要呈冲刷态势，河口河道呈淤积态势；与纯潮流作用下相比，河口区内48 h淤积变化量最高可达274.2%，淤积厚度高出一个量级.