人工采砂坑地形对水动力泥沙环境影响

谭海川，唐 亮，杨 洁，王 琮，张永强

（1. 深港天然气管道有限公司, 广东 深圳 518038；2. 中国石油集团工程技术研究有限公司, 天津 300451；3. 自然资源部 第一海洋研究所, 山东 青岛 266061）

海底管道是国家能源的主动脉，管道路由经常穿越不同的海洋功能区，容易受到海上生产作业的影响，人工采砂坑就是影响海底管道运行安全的重要风险点。深圳-香港海底管道支线（以下简称深港线）海底管道路由穿越人工海上采砂区，受到人工采砂、海流冲刷等因素的影响，采砂区不断扩大，对海底管道的在位稳定性造成严重影响。海底管道一旦发生损坏，影响不可估量。

针对海上人工采砂坑对海底管道路由区水动力泥沙环境的影响问题，国内外学者主要采用现场冲刷调查和数值模拟方法开展研究，取得了有益的结果。现场冲刷调查方法通常使用物探调查设备探测海底管道埋深和附近海床冲刷情况，对比分析调查结果，对管道悬空、裸露状态的变化和分布特征，找出演变规律和关键影响因素，从而为海底管道安全运行管理及悬跨段防护决策提供依据（吴明阳等, 2015; 金犇等, 2019; 沙红良等, 2023）。数值模拟方法是通过对现场观测数据的分析，建立海底管道路由一定范围内的水流泥沙数学模型，利用实测潮位、流速、流向及悬沙质量浓度资料率定并验证模型，进而模拟分析近年来人为开发活动对海底管道附近海床冲淤变化的影响，常用的数值模式包括DELFT3D、MIKE21 和FVCOM 等（曹祖德等, 2009; 英晓明等, 2014; 朱府等,2014; Rout et al, 2021; Song et al, 2022; Tseng et al, 2022）。

综上所述，现场实测资料的分析和数值模拟已经成为进行海床冲淤分析的有效方法。深港线海底管道路由区存在多处人工采砂坑，给管道在位稳定运行带来严重影响。本文针对深港线海底管道路由区KP5+458 点采砂坑形成后的冲淤问题，基于采砂坑及管道周边海域多期水深实测资料，采用实测资料对比和数值模拟方法，分析采砂坑存在对海底管道的潜在影响，以便制定有针对性的风险消减措施，保障深港线海底管道在位稳定运行。

1 研究区域

香港支线海底管道深圳段（深港线）起点为深圳大铲岛，终点为深港边界，海底段全长约14.5 km（图1）。采砂区位于海底管道西侧，整体呈盆状地形，水深为1～36 m，由东向西呈阶梯状分布。采砂区边坡顶部为原始浅滩，海床较平坦，水深约1～4 m；采砂坑边缘边坡陡峭，水深急剧增加到20 m，采砂区底部地势较为平坦，大部分水深在20 m 以深（图2）。

图1 研究区域位置Fig. 1 Location of the study area

图2 局部放大深港线采砂区状况Fig. 2 Local enlargement of the sand mining area along the Shenzhen-Hong Kong pipeline route

2018 年12 月，国家海洋局南海调查技术中心在对深港管道进行探测调查时发现，原有采砂区边缘有向管道路由方向扩展的趋势。2019 年3 月，南海调查中心进一步开展调查表明，该段管道西侧存在水深2～33 m 的采砂坑，形成平行管道路由长约6.7 km 的陡峭边坡，采砂坑距离海底管道最近处仅278 m，如图2 中KP5+458 点所示。采砂坑边坡坡度较大，人为或自然灾害易导致海床失稳和沙坑边坡滑塌，进而改变海管基础或导致管道裸露。目前采砂坑底部水深已经深于管道顶部标高，严重威胁管道安全。

2 研究方法

2.1 潮流场数值模型

采用MIKE 系列软件中的三角形网格水动力模块（FM 模块）计算潮流。FM 模块是软件核心的基础模块，可用于各种潮流场模拟和水动力学研究应用。其水流运动控制方程是二维浅水方程（DHI,2017a）：

式中：h为总水深，且h=η+d，其中 η和d分别为水面高度和静水深；t为时间；u¯ 和v¯分别为沿x和y方向的深度平均流速；S为点源流量；g为重力加速度；f为科氏参数； ρ为流体密度； ρ0为参考密度；us与vs为点源流速；Tij为应力项，包括黏性应力、紊流应力和对流等，根据水深平均的流速梯度计算。

在控制方程的求解过程中，FM 模块采用三角形网格，使用有限体积法进行离散；时间积分采用显式欧拉格式，计算速度较快，而且对于较小规模的潮流模拟问题，具有较高的计算效率；计算中采用干湿网格方法可以在模拟过程中动态调整边界条件和物理过程，能更准确地模拟潮汐的变化和浅滩的演化。

2.2 波浪场数值模型

采用MIKE 软件中SW 模块模拟波浪场，该模块适用于近岸海域、湖泊和河口的波浪数值模型，是国际上成熟、通用的海浪数值计算模式，包含空间的传播、因流场和水深场的改变而导致的波浪折射、波浪向浅滩的传播、波浪的反射和衍射、波浪生长和衰减过程、风浪的生长、波浪破碎对能量的消耗、波浪白帽的能量消耗、底摩擦效应，以及波-波相互作用等波浪传播过程。

在直角坐标系中，SW 模型控制方程可表示为（DHI, 2017b）：

式中：N为动谱密度；t为时间； ∇为微分算子； →v为波群速度；σ为相对频率；cx、cy、cθ和cσ分别为波群在四维相空间中的传播速度； θ为波向；s为沿 θ方向空间坐标；m为垂直于s的坐标； ∇→x为在x¯空间上的二维微分算子。

在分析采砂坑的地形冲淤变化和稳定性时，以年均波浪场作为动力场，首先须过滤掉对泥沙运动作用轻微的波浪（即波高小于0.5 m 的数据），统计各主导波向年平均波高，将其作为波要素进行波浪场反演，为泥沙数学模型提供波浪边界条件。

2.3 泥沙运动控制方程

采用MIKE 软件中的FM 模块模拟泥沙运动。FM 模块采用标准Galerkin 有限元法进行水平空间离散；在时间上，采用显式迎风差分格式离散动量方程与输运方程（DHI, 2017a）。泥沙控制方程为：

式中：c为沿深度平均的含沙量；u和v分别为沿x和y方向的深度平均流速；Dx和Dy分别为沿x和y方向的泥沙水平扩散系数，可取与水体紊动扩散系数相等的值；α为沉降几率或恢复饱和系数； ω为泥沙沉速，可根据平均底质粒径计算公式即张瑞瑾公式（周美蓉等, 2021）计算求得；C为床面冲淤强度；S*为波流共同作用下的挟沙量，根据波流挟沙的原理，S*可近似为：

式中S*C和S*W分别为潮流和波浪作用下的挟沙能量，可同时考虑潮流和波浪对泥沙的悬浮作用。

3 模型验证

为了解区域海区及采砂坑海域的动力环境、冲淤规律特征，采用波-流同步耦合模型建立目标海域潮流场和波浪场数值模型。在水动力模型基础上耦合泥沙冲淤模型，模拟目标海域潮流场、波浪场及泥沙冲淤环境，分析采砂坑及周边海域海底冲淤趋势。图3a 为广东近岸海域大区域模型范围及网格划分，对工程区网格进行加密处理；图3b 为珠江口伶仃洋和外部海域小区域模型范围；大区域数值模型采用二维/三维垂向平均模式，小区域数值模型采用三维垂向平均模式。

图3 模型计算范围及网格划分Fig. 3 Model calculation range and meshing

3.1 风场数据验证

研究区所在海域每年均会受台风过程影响，台风期间偏S 向风会导致内伶仃洋内大风增水和较大波浪传入。本文搜集赤湾、珠海、大万山及广州四个海洋站2017 年8 月的台风“天鸽”（编号1713），以及2018 年9 月的台风“山竹”（编号1822）期间的风速、风向实测资料，与模型采用的风场资料进行验证，具体站位分布见图1。对比分析实测数据与模拟数据（图4 和图5）可以看出，各海洋站实测风速和风向与模型采用的风速、风向资料基本吻合。

图4 海洋站计算风速和实测风速验证曲线Fig. 4 Verification curves of calculated and measured wind speed at ocean stations

图5 海洋站计算风向和实测风向验证曲线Fig. 5 Verification curves of calculated and measured wind direction at ocean stations

3.2 大区域潮流场模型潮位验证

本文搜集赤湾、珠海、大万山及广州四个海洋站2017 年8 月（1713 号台风“天鸽”）及2018 年9月（1822 号台风“山竹”）非台风影响期间的潮位实测资料，与模型计算潮位进行验证，结果表明，在各水文条件下，各潮位站模型与实测的潮位过程线吻合良好，模型的涨、落潮历时和相位与实测资料一致（图6）。

图6 海洋站计算潮位和实测潮位验证曲线Fig. 6 Verification curves of calculated and measured tide level at ocean station

3.3 小区域模型潮汐、潮流、波浪验证

利用研究海域C1～C5 五个潮流测站和C1、C5 两个临时潮位站大潮期的实测潮流资料对模型计算得到的潮流进行率定。C1、C5 潮位验证及C1～C5 测点流速流向验证结果分别见图7 和图8。

图7 大潮期潮位验证曲线Fig. 7 Verification curves of the tide level during the spring tide

从潮位验证结果来看，大、小潮各潮位站计算值与实测值基本一致。计算的高、低潮位出现的时间与实测高、低潮位出现的时间都吻合得较好。潮位平均偏差在0.04 m 左右，误差范围小于10%，满足JTS/T 231-2—2010《海岸与河口潮流泥沙模拟技术规程》（交通运输部水运局, 2010）的要求。从流速、流向验证结果来看，大、小潮各测点流速、流向计算值与实测值吻合较好，相位偏差较小，流速过程与现场基本一致，平均流速计算值与实测值偏差在10%左右，平均流向小于10°，也满足《海岸与河口潮流泥沙模拟技术规程》要求。因此，本模型计算得出的结果是合理的，模拟的流场反映了工程海域的潮流运动特征。

4 结果与分析

4.1 采砂区冲淤特征

自2018 年发现采砂坑以来，共开展了历时3 年（2019 至2021 年）的 8 次采砂坑及管线区水深地形测量和沉积物取样工作。现场调查资料充分，详细查清了采砂区的坑内淤积、边界后退、边坡变化情况，能够准确反映调查期间采砂坑冲淤发育趋势和规律。8 次水深地形调查结果如图9 所示。

为评估采砂坑的稳定性和地形演化情况，本文分析了KP4 和KP5 两个典型断面（断面位置见图1）自2019 年3 月到2022 年9 月的坡面角度变化及坡顶高程变化，结果如表1 和表2 所示。

表1 2019 年3 月至2022 年9 月坡面角度（°）变化Table 1 Changes in slope angle (°) from March 2019 to September 2022

表2 2019 年3 月至2022 年9 月坡顶高程（m）变化Table 2 Variation of slope crest elevation (m) from March 2019 to September 2022

对比图9 不同时期的水深地形可以看出：

1）2019 年3 月至2022 年9 月期间，采砂坑总体形态保持一致，未发生骤冲、骤淤现象。

2）采砂坑内总体为淤积环境，呈现淤平趋势，深坑淤积强度最大，其中北部深坑在2020 年4月前回淤明显，2021 年至今无明显冲淤现象；南部深坑总体为淤积现象，淤积变化不大。

3）2020 年9 月至2021 年1 月采砂坑北部深坑西侧、采砂坑中部发生局部变深，为采砂活动所致；2021 年1 月以后采砂坑内未发现采砂活动导致水深地形变化。

4）2020 年9 月前采砂区北部深坑淤积强度较大，与临近时间采砂活动关系密切。

整体看来，观测期间采砂坑形态总体保持不变，坑内淤积为主，无骤冲、骤淤现象，北部深坑2021 年1 月前强淤积过程与采砂活动关系密切，2021 年1 月以后采砂坑整体淤积变化不大。

图10 至图12 分别表示采砂区及周边海域1 a、2 a、5 a 后冲淤特征。通过对比分析，采砂坑的南北两侧及东侧为冲刷环境，局部冲刷强度为0.5 m/a 左右；坑内为淤积环境，局部最大淤积强度达到0.6 m/a 以上，整体淤积强度在0.2 m/a。长期条件下，采砂区的南北两侧及东侧为冲刷环境，局部冲刷深度分别约为1.1 m 和1.6 m，采砂区内为淤积环境，整体淤积强度分别约为0.5 m/a 和1 m/a。采砂区北侧边界向西北方向扩展，远离管道走向，对管道无影响。采砂区南侧边界向东南方向扩展，与管道走向有一定重叠，导致铜鼓水道附近KP12～KP14 段管道冲刷，局部强度分别为0.3 m/a、0.8 m/a左右，未达到管道埋设，其他管道段仍呈现淤积环境。

图10 采砂坑及周边海域1 a 后冲淤量分布Fig. 10 Distribution of scouring and silting amount in the sandpits and their surrounding waters after one year

图11 现状水深条件下2 a 后冲淤量分布Fig. 11 Distribution of scouring and silting amount under the case of present water depth after two years

图12 现状水深条件下5 a 后冲淤量分布Fig. 12 Distribution of scouring and silting amount under the case of present water depth after five years

4.2 采砂坑边界变化

模拟结果（图13）表明，采砂边界变化主要发生在采砂坑南北边界，东侧靠近管线段采砂坑边界变化略小。

图13 模拟采砂坑边界变化Fig. 13 Variations of the simulated boundary of the sand mining area

采砂坑南北边界后退主要受控于涨落潮海流垂直冲刷。1 a、2 a 和5 a 后，采砂区北侧边界分别向西北方向呈环状扩展0.4、0.8 和1.2 km，采砂区南侧边界向东南方向呈环状扩展0.5、0.9 和1.5 km。采砂区东侧边界后退主要由凸体局部冲刷所致，1 a、2 a 和5 a 后向靠近管道方向扩展6、11 和30 m，扩展段为与管道KP4+120 至KP11+120之间平行的部分。

4.3 采砂坑边坡变化

对于采砂区东侧边界的边坡，选取KP8+387、KP10+387 处2 条典型断面分析长期冲刷条件下的采砂区边坡变化特征，如图14 所示。可以看出，采砂坑边坡发育趋势为坡底淤积、边坡变缓、坡顶侵蚀。坡顶侵蚀会导致采砂坑边界线向管线扩展。

图14 长期冲刷作用下采砂区边坡变化Fig. 14 Slope changes of the sand mining area under the action of long-term scouring

5 结 论

本文分析了深港海底管道路由区域2019 年3 月至2022 年9 月期间8 次采砂坑及管线区水深地形测量资料，通过数值模拟并与实测资料相互验证，掌握区域海区及采砂坑海域的动力环境、冲淤规律特征。数值模拟方法能够准确反映调查期间采砂坑冲淤发育趋势和规律，分析结果可为海底管道路由采砂区制定有效防护方案、采取有效防护措施提供科学依据，其主要结论如下。

1）采砂坑总体形态保持不变，未发生骤冲、骤淤现象。

2）采砂坑内整体为淤积环境，采砂坑大坡度边坡主要分布在东南部边坡、南部边坡，以及东北侧局部；边坡坡度总体维持不变。

3）采砂坑东侧边缘走向与涨落潮主流向斜交，为冲刷环境，采砂坑边缘突出部分进行调整，形成凸体局部冲刷，采砂坑边缘线整体后退。

4）模拟冲刷1 年后东侧边界整体后退约6 m，模拟冲刷2 a 后东侧边界整体后退约11 m，模拟冲刷5 a 后东侧边缘整体后退约30 m。