大风江典型工程流场变化特征及悬浮泥沙扩散数值模拟*

陈逸航，高劲松

(南宁师范大学，北部湾环境演变与资源利用教育部重点实验室，广西南宁 530001)

大风江，又名平银江，上游(灵山伯劳河段)河面宽约40 m，平常水深0.8 m左右；
中游(那彭河段)河面宽约70 m，平常水深1.0 m左右；
下游(平银河段)河面宽约100 m，平常水深1.5 m左右。大风江流域有着天然的资源，不仅有国家级红树林保护区，还有丰富的渔业资源、浮游生物等。海洋工程施工期间所造成的悬浮泥沙扩散导致局部海域悬浮物增加，从而影响水体的透光性[1]，进而影响浮游植物的光合作用乃至海洋生态环境[2]。前人采用多种数值模拟方法分析了工程造成的悬浮泥沙扩散规律和特征[3-7]，为区域海洋生态环境保护提供了重要的科学依据和参考价值。针对北部湾悬浮泥沙研究，许晟轶[8]通过三维水动力数值模型计算了广西沿海重点港湾的悬浮泥沙分布。MIKE21模型为丹麦水力学研究所开发的二维水动力数学模型，广泛应用于海岸河口的潮流模拟和悬沙扩散输移[9-12]。一些学者通过MIKE21模型对广西北部湾典型港湾的悬浮泥沙扩散特征进行了详细分析[9,13-17]；
梁文等[18]运用遥感技术分析了廉州湾的悬浮泥沙动态特征。然而，针对大风江悬浮泥沙的研究报道较少见。罗亚飞等[19]通过实测数据结合Landsat TM、ETM卫星影像，对大风江附近海域悬沙进行遥感反演以及输移特征分析，发现大风江的悬沙浓度水平整体较小，江口以内水体悬沙浓度明显高于外海，且悬沙分布还受西南向沿岸流、季风等的影响。目前还未见有文章针对大风江的悬浮泥沙进行数值模拟分析。本研究以大风江大桥为例，基于MIKE模型分析流场变化特征以及悬浮泥沙扩散特征，为大风江的海洋环境保护提供科技支撑。

1.1 广西北部湾潮流模型构建与验证

为给大风江模型提供潮位边界条件，本研究先构建北部湾潮波数学模型。外海开边界由潮位控制，其潮位数据由NAO大洋潮汐预报模式提供，并结合沿岸潮位站验证资料予以调整。整体模式中，NAO99b及NAO99Jb提供了包含M2、S2、K1、O1、N2、P1、K2、Q1、M1、J1、OO1、2N2、Mu2、Nu2、L2、T2等共16个天文分潮的调和常数，适用于中国沿海特定期间的短期逐时潮位预报，控制方程组的数值求解采用有限体积法，模型计算区域采用非结构三角形网格，并在广西沿岸进行局部加密，网格间距局部岸线处约200 m，外海开边界最宽处约7 000 m，网格单元共计34 918个。

模型岸线采用现状岸线，越南一侧岸线采用美国海洋大气局(NOAA)提供的数据；
水深地形采用中国人民解放军海军司令部航海保证部2005年版之后的海图，广西沿岸局部港湾水深更新至2016年。水深及潮位资料统一至当地平均海平面。模型计算起止时间根据实测水文资料设定，时长约90 d。

为验证广西北部湾潮流模型的准确性，选取白龙尾、炮台角以及企沙潮位站2012年1月的潮位资料对模型进行验证，如图1所示(从上至下分别为白龙尾、炮台角以及企沙潮位站)，模拟和监测结果的相位变化和振幅都比较一致，验证结果较好，可为大风江潮流模型提供开边界。

图1 白龙尾、炮台角、企沙等3个观测站的监测水位和模拟水位验证Fig.1 Verification of monitoring water level and simulated water level at three observation stations of Bailongwei,Paotaijiao and Qisha

1.2 大风江二维潮流数值模型构建

大风江潮流数值模型的计算区域、工程区位置以及验证点布置如图2所示。计算范围为108.41°-109.17° E、21.39°-21.91° N，包括钦州湾、三娘湾和廉州湾，由于项目所在海域岛屿众多，岸线曲折，采用非结构三角形网格可以较好地贴合自然岸线，提高计算精度和计算效率，便于各种工程的准确布置。计算区域的工程前网格剖分如图2所示，在工程区附近进行网格加密。计算空间步长20-2 500 m，网格单元18 302个，网格节点10 004个。水深地形数据采用2012版流沙湾至东兴港海图、2012版钦州湾海图、2012版大风江口海图、广西908专项以及2019年8月局部海域调查数据，计算区域地形分布如图3所示。岸线采用最新的资料。外海潮位开边界由广西北部湾潮流模型提供，大风江上游开边界距工程区约16 km，本研究考虑径流量的作用。由于4月中旬为洪枯过渡季节，故大风江干流流量估取稍大于年平均径流量，即20 m3/s，丹竹江支流估取10 m3/s。水深及潮位均统一至当地平均海平面，坐标系统采用北京54坐标系。潮位验证资料时间为2020年4月14日11:00至2020年4月15日14:00，潮流实测资料为2020年4月14日12:00至2020年4月15日13:00在大风江海域大潮期调查的6个潮流站资料，时间与潮位站观测时间一致。模型计算时间从2020年4月1日至30日共30 d。

图2 大风江潮流模型计算区域及观测站分布图Fig.2 Dafengjiang River tidal current model calculation area and observation station distribution map

图3 计算区域网格划分图(a)和水深地形分布图(b)Fig.3 Grid partition map of calculation area (a) and distribution map of bathymetric topography (b)

计算拟建大桥桩基对流场的影响时，由于桩基的尺寸小于计算网格的尺度，网格无法反映桩基的存在，模拟计算参阅文献[20]对桥墩进行概化处理。

对桥墩所在网格的原水深进行水深修正：

式中，ΔH为桩基所在区的折减水深；
B1为桩墩迎水面宽度，B0为网格过水断宽度；
H为原网格水深。桥墩所在网格单元的糙率按下式计算局部阻力系数：

式中，n为桥墩所在单元的综合糙率；
nc为床面糙率；
nd为桥墩糙率；
ζ为桥墩阻力系数；
β为桥墩形状系数，取0.85。

1.3 大风江二维潮流数值模型验证

对三娘湾潮位站实测水位过程与计算值进行比较(图4)，起始时间为2020年4月14日11：00，图4中红线为计算的潮位值，黑点“◆”为实测值。从图4可以看出，计算的潮位过程与实测资料吻合较好。验证结果表明采用的二维潮流数值模型能较好地模拟大风江所在海域水位变化过程，也为准确模拟当地的潮流变化过程奠定基础。

图4 三娘湾潮位站的监测水位与模拟水位验证Fig.4 Verification of monitoring water level and simulation water level at the Sanniang Bay tide station

图5给出了4个潮流测站(1#、2#、3#、4#)的流速计算结果与实测结果比较，流向以北方向为起始，顺时针旋转为正。由图5可知，各验证点计算流速和实测资料基本吻合，流向验证较好，最大误差小于10%，总体来看，流速过程线的形态基本一致，表明建立的二维潮流数学模型能较好地模拟工程所在海区水流传播过程和水流运动规律。

图5 4个站点的流速和流向验证Fig.5 Flow velocity and direction verification at 4 stations

1.4 悬浮泥沙扩散模型

1.4.1 基本方程

采用二维悬沙输运方程预测施工期间产生的悬浮物对水质的影响，平面二维悬沙运动方程如下：

式中，C为垂向平均含沙量，ε为垂向平均的扩散系数，

FC=SC+

式中，SC为输入源强，a为沉积系数，M为冲刷系数，τb为底部切应力，τe为临界冲刷切应力，τd为临界淤积切应力。通过联立水动力方程数值求解悬浮物扩散方程。

1.4.2 源强计算

废水泥沙的产生量与管桩下压的深度、管桩体积和施工抽水工况等因素有关，其进入海洋环境的泄漏量可按产生量的5%估算。钢护筒内废水泥沙的产生量采取以下公式进行测算：

G=1/4πd2·h·ρ，

其中,G为桩基施工时产生的护筒内泥沙量，分为不同桥段相应泥沙产生量，此处取大桥区段平均值;d为护筒直径，比桩基本身略大10-20 cm，本研究取15 cm，本工程桥墩桩基直径为2.20 m，引桥墩1.80 m，则d桥墩桩为2.35 m，d引桥墩桩为1.95 m；
h为各区段海底覆盖层厚度，平均约10 m；
ρ为覆盖层泥沙浓度，约为1.47×103kg/m3。

根据上述公式，单日(施工10 h计)悬浮物泄漏量按照泥沙产生量的5%估算，则本工程单个主桥墩桩基施工产生悬浮物泄漏源单点强为0.088 6 kg/s；
每组主桥墩同一时间只能施工2个桩基，则悬浮物源强为0.177 2 kg/s；
每组引桥墩同一时间只能施工1个桩基，则单个引桥墩桩基施工产生悬浮物泄漏源单点强为0.061 0 kg/s。本研究选择低潮时作为泥沙的释放时间。

2.1 大风江潮流场特征

由图6和图7可知，浅滩、岛屿周围以及岸边流速相对较小，流向多变；
航道和深槽处流速较大，流向与航道、深槽走向基本一致。涨急时外海水沿着西北向进入大风江，并顺着河道上溯至上游地区，最大流速为0.80 m/s；
落急时，上游各支汊水流汇入主流后沿东南偏南向到达河口地区，再转为西南向流向外海，最大流速为0.89 m/s。大风江海域落急流速(平均流速为0.42 m/s)大于涨急流速(平均流速为0.28 m/s)。拟建的大桥区域潮流流向受地形限制，涨急和落急方向为西北-东南向。

图6 工程前大风江海域涨急和落急时刻流场Fig.6 Flow field of flood and ebb peak in the Dafengjiang River sea area before construction

图7 工程前工程区附近局部海域涨急和落急流场Fig.7 Flow field of flood and ebb peak in local sea area near project area before construction before construction

2.2 大桥建设对大风江潮流场影响

本研究模拟了大桥建设后的涨急和落急流场分布(图8)，发现整个大风江的潮流场变化不大，流场的变化主要集中在大桥附近。为定量分析大桥建设对大风江的潮流场影响情况，在拟建大桥周边海域选取20个特征点(T1-T20，图9)。由表1可知，无论是涨急还是落急时刻，大桥建设对流场的影响主要集中在桥两端2 km距离范围内。涨急时刻，特征点T5和T12的流速变化在6%以上；
落急时刻，特征点T12和T13的流速变化在8%以上，集中在桥的东北端。桥的西南端(T1、T8、T9、T16和T17)流速变化较小。因此，桥墩对流场的影响主要在桥的东北端海域。

图8 大桥建设后附近海域涨急和落急流场Fig.8 Flow field of flood and ebb peak near sea area after bridge construction

图9 大桥附近特征点分布Fig.9 Distribution of feature points around the bridge

表1 大桥建设前后特征点流速比较Table 1 Comparison of flow velocity of feature points before and after bridge construction

续表Continued table

2.3 悬浮泥沙扩散特征

基于潮流模型，本研究对拟建大桥桥墩施工引起的悬浮泥沙扩散特征进行模拟(图10)，统计结果如表2所示，泥沙释放时间为低潮时。由于落潮流大于涨潮流，因此落潮期间浓度10 mg/L的悬浮泥沙的包络面积(10.11 km2)远大于涨潮期间的包络面积(1.10 km2)。落潮时最远扩散距离为7.15 km，涨潮时最远扩散距离为2.69 km。悬浮泥沙的扩散方向与涨潮流和落潮流方向一致，均沿着西北向和东南向扩散。由于悬浮泥沙的扩散会导致局部海域悬浮物增加，从而影响水体的透光性[1]，进而影响浮游植物的光合作用乃至海洋生态环境[2]，因此建议控制施工规模，并采取防污帘等措施来减少悬浮泥沙的扩散。

图10 涨潮和落潮期间悬浮泥沙浓度>10 mg/L的包络面积Fig.10 Envelope area of suspended sediment with the concemtration >10 mg/L during flood and ebb tide

表2 涨潮和落潮期间悬浮泥沙增量包络面积及最远扩散距离Table 2 Incremental envelope area and the largest diffusion distance of suspended sediment during flood tide and ebb tide

海洋工程施工期间所造成的悬浮泥沙扩散导致局部海域悬浮物增加，从而影响水体的透光性及海洋生态环境。本研究通过构建广西北部湾潮流模型从而提供开边界数据，进而构建大风江潮流模型和悬浮泥沙扩散模型，发现模拟的潮位和潮流结果与实测数据吻合较好。大风江呈现往复流特征，涨急时最大流速为0.80 m/s，落急时最大流速为0.89 m/s，落急流速(平均流速为0.42 m/s)大于涨急流速(平均流速为0.28 m/s)。拟建大桥建设对大风江的潮流场影响有限，影响主要集中在桥的东北端海域。当悬浮泥沙在低潮释放时，落潮期间浓度>10 mg/L的悬浮泥沙包络面积(10.11 km2)远大于涨潮期间的包络面积(1.10 km2)。同时落潮时最远扩散距离为7.15 km，涨潮时最远扩散距离为2.69 km。悬浮泥沙的扩散方向与涨潮流和落潮流方向一致，分别沿着西北向和东南向扩散。悬浮泥沙的扩散会对周边海域的水质造成一定影响，因此建议控制施工规模并采取防污帘等措施来减少悬浮泥沙扩散。