黄河三角洲海岸冲淤及泥沙输运模型数值模拟分析

王奎峰，张太平，宋新强，许国辉,尚桂勇,郑建国

(1.山东省地质科学研究院，国土资源部金矿成矿过程与资源利用重点实验室，山东省金属矿产成矿地质过程与资源利用重点实验室，山东 济南 250013；2.山东大学土建与水利学院，山东 济南 250061；3.山东省沉积成矿作用与沉积矿产重点实验室，山东 青岛 266590；4.东营市国土资源局河口分局，山东 东营 257100；5.中国海洋大学环境科学与工程学院，山东 青岛 266100；6.东营市公路管理局河口分局，山东 东营 257200)

河口是河流沉积物向海域传输的通道，河口泥沙传输过程相当复杂，河口泥沙输运过程对于研究河口陆海相互作用及近岸地质环境演化有重要的参考意义[1-2]。作为最新成陆的黄河三角洲地区，陆海相互作用强烈，1855 年以来，形成了以宁海为顶点的扇形现代黄河三角洲[3-8]。超过80%的黄河入海泥沙沉积在河口三角洲区域，而河口动力过程则是形成这种三角洲沉积格局的重要原因。黄河携巨量泥沙入海并在河口地区快速落淤，对黄河口海域水动力环境和近海生态环境影响深远[1，9]。因此，研究黄河口的泥沙运动、冲淤动态及其效应对于黄河口治理及可持续发展具有重要意义。以往人们对于河口海岸区域海洋动力和悬浮泥沙输运的研究主要依赖于观测，但现场观测的经济成本较高，同时也不能做到时空上的连续，且存在极端天气条件不利影响因素，效果不好[1，10-14]。该文结合近年来在该地区开展的基金及地勘项目成果数据，采用三维水动力—泥沙输运模型对黄河三角洲区域海岸冲淤和悬浮泥沙输送进行数值模拟，对认识泥沙输运、盐度分布，河口混合过程等具有重要意义[5,15-17]。

1 研究区概况

黄河三角洲地区是具有国家战略地位的高效经济区，是一典型扇形三角洲，由黄河填海造陆而形成，属河流冲积物覆盖海相层的二元相结构，为世界上最年轻的陆地之一[18]。黄河口区域河流动力与海洋动力相互作用强烈，入海水沙变化幅度较大，泥沙在多种动力因子控制下的运动机制相当复杂。近年来，由于人类工程经济活动及自然因素的影响，黄河入海水沙量急剧减小[1，19]。同时，黄河三角洲海岸带地形地貌及岸线的变化也使得黄河口海域水动力条件及泥沙输运发生了变化。

图1 黄河三角洲区域位置图

2 数学模型简介

河口海岸数学模型是研究河口沉积动力过程的有效手段之一。该文采用的三维河口海岸动力模型HEM-3D是由美国佛吉尼亚海洋研究所(VIMS)开发的，是一个用于研究河流、湖泊、水库、河口及陆架海区的水动力、物质输运和生物地球化学过程的数学模型包[16]。HEM-3D模型已经在国内外多个研究项目中得到了广泛的应用，效果良好，非常适用于黄河口多环境动力因子影响下河口泥沙输运的研究。

2.1 模型基本结构

该文主要针对HEM-3D模型水动力模块和泥沙输运模块开展应用研究。HEM-3D的水动力模型考虑河口海岸区域淡水径流、气象作用力、地形和底摩擦作用等输入因素影响下的潮流动力场，求解水动力场的自由表面高程、流速、湍流混合、盐度和温度等要素的时空变化。HEM-3D模型的泥沙输运在海底边界上考虑泥沙的沉降和再悬浮过程及悬沙浓度对泥沙沉速的影响，在水平边界上考虑随时间变化的物质输入和输出过程方程忽略水平的湍扩散项，采用和盐度、温度输运相同的高阶对流—扩散格式求解[2]。

2.2 数学模型设置

(1)网格配置

为了准确拟合模型区域的几何边界，数值模拟采用曲线—正交网格来较为准确地处理模型区域的边界，岸线边界区采用三角形网格。模型在垂向上采用Sigma坐标变换，分两层进行计算。

(2)水深

研究区域的水深资料主要采用2013年36条实测断面的水深数据(图2)。结合收集历史水深资料，经过校准和率定后，插值到各个计算网格的中心点上，为数学模型提供水深条件。这些资料数据主要来源为黄河口水文水资源勘测局在该区域内常年布设的监测断面数据。

图2 黄河三角洲实测断面位置图

(3)开边界设置

①外海开边界条件：在开边界上采用M2，S2，K1,O1四个主要分潮的调和常数来确定水位的变化。

式中：ζn，τn，Tn分别为每个分潮的振幅，初相位和周期；ζser为余水位的时间序列。开边界上各分潮的调和常数由渤海潮汐数学模型计算结果插值给出(图2)。

图3 渤海海域M2分潮同潮图

②海面边界条件：在海面z=1处

运动学边界条件：w|z=1=0

动力学边界条件：

海面的边界条件主要考虑风应力作用。其中，海面风应力项τsx和τsv由下式确定：

式中：Uw和Vw分别为海面10m高度上的风速在曲面-正交坐标x，y方向上的分量，风应力系数由下式确定：

式中：ρa和ρw分别为空气和水体的密度。

③海底边界条件：海底z=0处的运动学边界条件为垂向速度为零，w|z=0。

动力学边界条件由底部剪切应力控制：

底部相对粗糙高度z0，根据Heathersaw的研究，对于淤泥质海区，底部绝对粗糙高度在(2～7)×10-4m的范围内。

④海岸边界条件：采用滑动边界条件，边界处的法向流速为零，盐度、泥沙不存在通量，即：

开边界上的悬浮泥沙根据历史观测资料给定一个相对稳定的浓度值，孤东海区的悬浮泥沙主要考虑为潮流引起的再悬浮，现行黄河口泥沙经过再悬浮和潮流输运也可以抵达该区，北侧钓口三角洲区域风浪再悬浮的泥沙也对该海区的悬浮泥沙有主要影响。泥沙沉降速度的确定采用经验关系[2]：

Ws=Cαexp(-4.21+0.147G)

3 数值模拟计算结果分析

采用HEM-3D三维水动力和泥沙输运数学模型，模拟了黄河口海域的潮流场、盐度场、悬沙浓度场和近海域冲淤动态分布特征。

3.1 潮流场

该海区的潮流基本上属于规则半日潮流，为沿岸的往复流，涨潮流由西北向南，落潮流由南向西北，流速大致为20～80cm/s左右，莱州湾区域和渤海湾区域的涨落潮流存在大约6h的时间差。对模式输出的潮流时间序列(t=0～t=5/6T)进行分析，计算得到一个潮周期平均潮流场(图4)。涨潮时，潮流由西北向南，在河口的南侧形成回旋流，流速减小。黄河三角洲沿岸区域存在2个高流速中心，最大流速约为1.1m/s，分别位于三角洲北部(废弃刁口三角洲叶瓣前缘)和现行河口口门前缘[19]。

根据孤东区域的有关观测数据，选取GD01点位的潮流资料与模型结果对比，模拟结果显示，孤东海区的潮流基本上为东南(涨潮)—西北(落潮)的沿岸往复流性质，转流时间很短，且落潮流速高于涨潮流，该区域修建的沿岸人工海堤导致海岸线走向多变，形成潮流场与海岸相互作用强烈，海底剪切应力加强，模型计算结果与实测潮流基本吻合(图5)。

图4 t=0～t=5/6T时刻海域潮流场图

图5 实测潮流与模型结果对比图(GD01)

图6 t=0～t=5/6T时刻海域盐度分布图

图7 t=0～t=5/6T时刻海域悬浮泥沙浓度分布图

3.2 盐度分布

海域的盐度分布特征整体表现为三角洲的北部盐度较高，在18～30PSU之间，孤东以南的海域盐度较低，尤其是现行河口口门附近区域，由于受河流入海冲淡水的影响，盐度最低，莱州湾内的盐度大致在15PSU左右，在远离河口的区域盐度升高。另外，三角洲的近岸区域和外海区域的盐度也有差别。在潮周期内，三角洲现行河口区域的盐度随潮流场发生变化，而孤东以北的区域相对稳定。现行河口口门区域是径流和潮流相互作用强烈的区域，河口羽状流的向海扩展受控于涨落潮流的流速大小和潮流方向，导致羽状流的路径和向海扩展的幅度随时间变化(图6)。

3.3 悬浮泥沙输运

三角洲海域的悬浮泥沙浓度分布与潮流场变化和河口泥沙输入有密切的关系(图7)。在空间上存在多源控制的特征。在三角洲北部(刁口三角洲叶瓣前缘)，受五号桩外强潮流区的影响，近岸海底的泥沙发生明显的再悬浮，并在涨潮流向南输送，含沙量达1.5g/L左右。在三角洲南部(现行河口区域和莱州湾区域)，受现行河口入海泥沙扩散的影响显著。现行河口的口门区域可以看到明显的河流泥沙随羽状流向海扩展的分布形态，从河口入海后向南偏转，悬沙浓度约为3.5g/L。1976—1996年行河期间形成的河口叶瓣，由于地形向海突出，与潮流场相互作用强烈，海底泥沙再悬浮明显，在莱州湾北侧形成一个明显的高泥沙浓度中心，而莱州湾的近岸区域泥沙浓度较低。数值模拟显示的三角洲海域悬浮泥沙浓度的空间分布形态与卫星遥感数据基本吻合(图8)。

从潮周期内的变化来看，黄河三角洲区域始终存在再悬浮，泥沙浓度近岸高于深水区域；现行河口的羽状流随潮流变化而变化，羽状流的路径以及扩展幅度变化，引起了盐度和悬浮泥沙浓度均发生调整。莱州湾北侧的高浓度中心在潮周期内始终存在。

图8 2013年9月11日LANDSAT卫星显示黄河三角洲悬浮泥沙浓度高值区图

3.4 近海岸冲淤分布

根据黄河三角洲海域年冲淤分布的数值模拟计算结果(图9)，在北部废弃三角洲刁口三角洲叶瓣区域，由于海洋动力作用强烈，海底泥沙再悬浮活跃，再悬浮泥沙在斜坡重力牵引下向深水区输运，导致该区域浅水冲刷，海底侵蚀显著，形成了明显呈沿岸展布的侵蚀中心，中心的最大侵蚀厚度约为0.25m左右，位于离岸5～7km左右，在侵蚀中心以外，侵蚀快速减弱，在离岸约20km以外形成弱微淤积，淤积厚度不足0.1m，淤积带有向渤海中部延伸的趋势，可能对渤海中部的泥质沉积物有密切的关系。在现行河口区域，以淤积为主，河口口门的最大淤积厚度约为0.4m，主要取决于黄河入海的泥沙量以及在海洋动力作用下的扩散过程。在羽状流扩散的控制下，泥沙入海后向南输运，因此河口口门的淤积中心向南延伸，至莱州湾区域逐渐减弱，外缘的淤积厚度为0.1m左右。海域的其他区域冲淤不明显[6，14,19-20]。这与李殿魁牵头联合黄河水利委员会黄河河口研究院及中科院地理所等单位开展的“黄河三角洲国土防护与生态修复技术研究”的研究成果[21]基本一致，也验证了数值模拟计算结果可行性[21]。

图9 黄河三角洲海域年冲淤量分布图(m)

与以往的研究文献不同的是，该次还选取了黄河三角洲的4个关键断面(图10)进行冲淤分析：A-A断面位于河口区新户镇北部冲刷中心的外缘，在废弃刁口三角洲叶瓣的西侧，在近岸15km以内，以冲刷为主，最大冲刷厚度不足0.1m，10km以外为淤积，随着水深加大，淤积厚度加大，深水区的淤积厚度约为0.06m(图11)；D-D断面位于东营市临港产业区，在黄河海港的东北侧，处于废弃刁口三角洲叶瓣冲刷带，近岸10km以内为冲刷，年最大的冲刷深度为0.2m，15km以外形成淤积，主要是由于近岸再悬浮泥沙在斜坡重力的牵引作用下向深水区输运并形成淤积，但是年淤积厚度不足0.1m(图12)；E-E断面位于黄河海港南侧，为废弃刁口三角洲叶瓣冲刷区的外侧，近岸8km以内为冲刷，年最大冲刷厚度约为0.11m，在8km以外转为淤积，年淤积厚度较小，不足0.05m(图13)；F-F断面位于莱州湾东侧的羊口镇，远离现行河口的淤积中心，同时莱州湾内动力较弱，泥沙运动不活跃，整体处于冲淤平衡状态，从冲淤分布看，近岸5km为冲刷，年冲刷厚度不足0.05m，外侧为淤积或冲淤平衡状态，趋势不明显(图14)，这与东营市黄河口泥沙研究所开展的监测断面水下地形冲淤实测数据结果[19-23]也大致相同[注]王奎峰．黄河三角洲高效生态经济区人工海岸生态地质环境调查报告，2016年。，验证了该次模拟分析的适用性和先导性。

图10 黄河三角洲关键断面位置图

图11 A-A断面的冲淤分布图

图13 E-E断面的冲淤分布图

图14 F-F断面的冲淤分布图

4 结论

该文运用HEM-3D模型，采用最新的黄河三角洲断面监测数据，针对黄河口三角洲的流场和泥沙扩散进行了三维数值模拟，计算结果基本反映了黄河口水动力场、悬浮泥沙、盐度和冲淤的变化特征。

(1)黄河三角洲潮流基本上属于规则半日潮流，为沿岸的往复流，涨潮流由西北向南，落潮流由南向西北，黄河三角洲沿岸区域存在2个高流速中心，分别位于三角洲北部和现行河口口门前缘。

(2)海域的盐度分布特征整体表现为三角洲的北部盐度较高，在18～30PSU之间，孤东以南的海域盐度较低，现行河口口门附近盐度最低，莱州湾内的盐度大致在15PSU左右，在远离河口的区域盐度升高。另外，三角洲的近岸区域和外海区域的盐度也有差别。在潮周期内三角洲现行河口区域的盐度随潮流场发生变化，而孤东以北的区域相对稳定。

(3)三角洲海域的悬浮泥沙浓度分布与潮流场变化和河口泥沙输入有密切的关系。在三角洲北部受五号桩外强潮流区的影响，近岸海底的泥沙发生明显的再悬浮，并在涨潮流向南输送，含沙量达1.5g/L左右。在三角洲南部(现行河口区域和莱州湾区域)，受现行河口入海泥沙扩散的影响显著。

(4)近海域冲淤分布方面，黄河三角洲北部刁口区域海底冲刷侵蚀严重，形成了明显呈沿岸展布的侵蚀中心，向南侵蚀逐渐减弱。在现行河口区域，以淤积为主，河口口门的最大淤积厚度约为0.4m，河口口门的淤积中心向南延伸，至莱州湾区域逐渐减弱。