人工采砂坑地形對(duì)水動(dòng)力泥沙環(huán)境影響

譚海川，唐 亮，楊 潔，王 琮，張永強(qiáng)

（1. 深港天然氣管道有限公司, 廣東 深圳 518038；2. 中國(guó)石油集團(tuán)工程技術(shù)研究有限公司, 天津 300451；3. 自然資源部 第一海洋研究所, 山東 青島 266061）

海底管道是國(guó)家能源的主動(dòng)脈，管道路由經(jīng)常穿越不同的海洋功能區(qū)，容易受到海上生產(chǎn)作業(yè)的影響，人工采砂坑就是影響海底管道運(yùn)行安全的重要風(fēng)險(xiǎn)點(diǎn)。深圳-香港海底管道支線（以下簡(jiǎn)稱(chēng)深港線）海底管道路由穿越人工海上采砂區(qū)，受到人工采砂、海流沖刷等因素的影響，采砂區(qū)不斷擴(kuò)大，對(duì)海底管道的在位穩(wěn)定性造成嚴(yán)重影響。海底管道一旦發(fā)生損壞，影響不可估量。

針對(duì)海上人工采砂坑對(duì)海底管道路由區(qū)水動(dòng)力泥沙環(huán)境的影響問(wèn)題，國(guó)內(nèi)外學(xué)者主要采用現(xiàn)場(chǎng)沖刷調(diào)查和數(shù)值模擬方法開(kāi)展研究，取得了有益的結(jié)果?，F(xiàn)場(chǎng)沖刷調(diào)查方法通常使用物探調(diào)查設(shè)備探測(cè)海底管道埋深和附近海床沖刷情況，對(duì)比分析調(diào)查結(jié)果，對(duì)管道懸空、裸露狀態(tài)的變化和分布特征，找出演變規(guī)律和關(guān)鍵影響因素，從而為海底管道安全運(yùn)行管理及懸跨段防護(hù)決策提供依據(jù)（吳明陽(yáng)等, 2015; 金犇等, 2019; 沙紅良等, 2023）。數(shù)值模擬方法是通過(guò)對(duì)現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)的分析，建立海底管道路由一定范圍內(nèi)的水流泥沙數(shù)學(xué)模型，利用實(shí)測(cè)潮位、流速、流向及懸沙質(zhì)量濃度資料率定并驗(yàn)證模型，進(jìn)而模擬分析近年來(lái)人為開(kāi)發(fā)活動(dòng)對(duì)海底管道附近海床沖淤變化的影響，常用的數(shù)值模式包括DELFT3D、MIKE21 和FVCOM 等（曹祖德等, 2009; 英曉明等, 2014; 朱府等,2014; Rout et al, 2021; Song et al, 2022; Tseng et al, 2022）。

綜上所述，現(xiàn)場(chǎng)實(shí)測(cè)資料的分析和數(shù)值模擬已經(jīng)成為進(jìn)行海床沖淤分析的有效方法。深港線海底管道路由區(qū)存在多處人工采砂坑，給管道在位穩(wěn)定運(yùn)行帶來(lái)嚴(yán)重影響。本文針對(duì)深港線海底管道路由區(qū)KP5+458 點(diǎn)采砂坑形成后的沖淤問(wèn)題，基于采砂坑及管道周邊海域多期水深實(shí)測(cè)資料，采用實(shí)測(cè)資料對(duì)比和數(shù)值模擬方法，分析采砂坑存在對(duì)海底管道的潛在影響，以便制定有針對(duì)性的風(fēng)險(xiǎn)消減措施，保障深港線海底管道在位穩(wěn)定運(yùn)行。

1 研究區(qū)域

香港支線海底管道深圳段（深港線）起點(diǎn)為深圳大鏟島，終點(diǎn)為深港邊界，海底段全長(zhǎng)約14.5 km（圖1）。采砂區(qū)位于海底管道西側(cè)，整體呈盆狀地形，水深為1～36 m，由東向西呈階梯狀分布。采砂區(qū)邊坡頂部為原始淺灘，海床較平坦，水深約1～4 m；采砂坑邊緣邊坡陡峭，水深急劇增加到20 m，采砂區(qū)底部地勢(shì)較為平坦，大部分水深在20 m 以深（圖2）。

圖1 研究區(qū)域位置Fig. 1 Location of the study area

圖2 局部放大深港線采砂區(qū)狀況Fig. 2 Local enlargement of the sand mining area along the Shenzhen-Hong Kong pipeline route

2018 年12 月，國(guó)家海洋局南海調(diào)查技術(shù)中心在對(duì)深港管道進(jìn)行探測(cè)調(diào)查時(shí)發(fā)現(xiàn)，原有采砂區(qū)邊緣有向管道路由方向擴(kuò)展的趨勢(shì)。2019 年3 月，南海調(diào)查中心進(jìn)一步開(kāi)展調(diào)查表明，該段管道西側(cè)存在水深2～33 m 的采砂坑，形成平行管道路由長(zhǎng)約6.7 km 的陡峭邊坡，采砂坑距離海底管道最近處僅278 m，如圖2 中KP5+458 點(diǎn)所示。采砂坑邊坡坡度較大，人為或自然災(zāi)害易導(dǎo)致海床失穩(wěn)和沙坑邊坡滑塌，進(jìn)而改變海管基礎(chǔ)或?qū)е鹿艿缆懵?。目前采砂坑底部水深已?jīng)深于管道頂部標(biāo)高，嚴(yán)重威脅管道安全。

2 研究方法

2.1 潮流場(chǎng)數(shù)值模型

采用MIKE 系列軟件中的三角形網(wǎng)格水動(dòng)力模塊（FM 模塊）計(jì)算潮流。FM 模塊是軟件核心的基礎(chǔ)模塊，可用于各種潮流場(chǎng)模擬和水動(dòng)力學(xué)研究應(yīng)用。其水流運(yùn)動(dòng)控制方程是二維淺水方程（DHI,2017a）：

式中：h為總水深，且h=η+d，其中 η和d分別為水面高度和靜水深；t為時(shí)間；uˉ 和vˉ分別為沿x和y方向的深度平均流速；S為點(diǎn)源流量；g為重力加速度；f為科氏參數(shù)； ρ為流體密度； ρ0為參考密度；us與vs為點(diǎn)源流速；Tij為應(yīng)力項(xiàng)，包括黏性應(yīng)力、紊流應(yīng)力和對(duì)流等，根據(jù)水深平均的流速梯度計(jì)算。

在控制方程的求解過(guò)程中，F(xiàn)M 模塊采用三角形網(wǎng)格，使用有限體積法進(jìn)行離散；時(shí)間積分采用顯式歐拉格式，計(jì)算速度較快，而且對(duì)于較小規(guī)模的潮流模擬問(wèn)題，具有較高的計(jì)算效率；計(jì)算中采用干濕網(wǎng)格方法可以在模擬過(guò)程中動(dòng)態(tài)調(diào)整邊界條件和物理過(guò)程，能更準(zhǔn)確地模擬潮汐的變化和淺灘的演化。

2.2 波浪場(chǎng)數(shù)值模型

采用MIKE 軟件中SW 模塊模擬波浪場(chǎng)，該模塊適用于近岸海域、湖泊和河口的波浪數(shù)值模型，是國(guó)際上成熟、通用的海浪數(shù)值計(jì)算模式，包含空間的傳播、因流場(chǎng)和水深場(chǎng)的改變而導(dǎo)致的波浪折射、波浪向淺灘的傳播、波浪的反射和衍射、波浪生長(zhǎng)和衰減過(guò)程、風(fēng)浪的生長(zhǎng)、波浪破碎對(duì)能量的消耗、波浪白帽的能量消耗、底摩擦效應(yīng)，以及波-波相互作用等波浪傳播過(guò)程。

在直角坐標(biāo)系中，SW 模型控制方程可表示為（DHI, 2017b）：

式中：N為動(dòng)譜密度；t為時(shí)間； ?為微分算子； →v為波群速度；σ為相對(duì)頻率；cx、cy、cθ和cσ分別為波群在四維相空間中的傳播速度； θ為波向；s為沿 θ方向空間坐標(biāo)；m為垂直于s的坐標(biāo)； ?→x為在xˉ空間上的二維微分算子。

在分析采砂坑的地形沖淤變化和穩(wěn)定性時(shí)，以年均波浪場(chǎng)作為動(dòng)力場(chǎng)，首先須過(guò)濾掉對(duì)泥沙運(yùn)動(dòng)作用輕微的波浪（即波高小于0.5 m 的數(shù)據(jù)），統(tǒng)計(jì)各主導(dǎo)波向年平均波高，將其作為波要素進(jìn)行波浪場(chǎng)反演，為泥沙數(shù)學(xué)模型提供波浪邊界條件。

2.3 泥沙運(yùn)動(dòng)控制方程

采用MIKE 軟件中的FM 模塊模擬泥沙運(yùn)動(dòng)。FM 模塊采用標(biāo)準(zhǔn)Galerkin 有限元法進(jìn)行水平空間離散；在時(shí)間上，采用顯式迎風(fēng)差分格式離散動(dòng)量方程與輸運(yùn)方程（DHI, 2017a）。泥沙控制方程為：

式中：c為沿深度平均的含沙量；u和v分別為沿x和y方向的深度平均流速；Dx和Dy分別為沿x和y方向的泥沙水平擴(kuò)散系數(shù)，可取與水體紊動(dòng)擴(kuò)散系數(shù)相等的值；α為沉降幾率或恢復(fù)飽和系數(shù)； ω為泥沙沉速，可根據(jù)平均底質(zhì)粒徑計(jì)算公式即張瑞瑾公式（周美蓉等, 2021）計(jì)算求得；C為床面沖淤?gòu)?qiáng)度；S*為波流共同作用下的挾沙量，根據(jù)波流挾沙的原理，S*可近似為：

式中S*C和S*W分別為潮流和波浪作用下的挾沙能量，可同時(shí)考慮潮流和波浪對(duì)泥沙的懸浮作用。

3 模型驗(yàn)證

為了解區(qū)域海區(qū)及采砂坑海域的動(dòng)力環(huán)境、沖淤規(guī)律特征，采用波-流同步耦合模型建立目標(biāo)海域潮流場(chǎng)和波浪場(chǎng)數(shù)值模型。在水動(dòng)力模型基礎(chǔ)上耦合泥沙沖淤模型，模擬目標(biāo)海域潮流場(chǎng)、波浪場(chǎng)及泥沙沖淤環(huán)境，分析采砂坑及周邊海域海底沖淤趨勢(shì)。圖3a 為廣東近岸海域大區(qū)域模型范圍及網(wǎng)格劃分，對(duì)工程區(qū)網(wǎng)格進(jìn)行加密處理；圖3b 為珠江口伶仃洋和外部海域小區(qū)域模型范圍；大區(qū)域數(shù)值模型采用二維/三維垂向平均模式，小區(qū)域數(shù)值模型采用三維垂向平均模式。

圖3 模型計(jì)算范圍及網(wǎng)格劃分Fig. 3 Model calculation range and meshing

3.1 風(fēng)場(chǎng)數(shù)據(jù)驗(yàn)證

研究區(qū)所在海域每年均會(huì)受臺(tái)風(fēng)過(guò)程影響，臺(tái)風(fēng)期間偏S 向風(fēng)會(huì)導(dǎo)致內(nèi)伶仃洋內(nèi)大風(fēng)增水和較大波浪傳入。本文搜集赤灣、珠海、大萬(wàn)山及廣州四個(gè)海洋站2017 年8 月的臺(tái)風(fēng)“天鴿”（編號(hào)1713），以及2018 年9 月的臺(tái)風(fēng)“山竹”（編號(hào)1822）期間的風(fēng)速、風(fēng)向?qū)崪y(cè)資料，與模型采用的風(fēng)場(chǎng)資料進(jìn)行驗(yàn)證，具體站位分布見(jiàn)圖1。對(duì)比分析實(shí)測(cè)數(shù)據(jù)與模擬數(shù)據(jù)（圖4 和圖5）可以看出，各海洋站實(shí)測(cè)風(fēng)速和風(fēng)向與模型采用的風(fēng)速、風(fēng)向資料基本吻合。

圖4 海洋站計(jì)算風(fēng)速和實(shí)測(cè)風(fēng)速驗(yàn)證曲線Fig. 4 Verification curves of calculated and measured wind speed at ocean stations

圖5 海洋站計(jì)算風(fēng)向和實(shí)測(cè)風(fēng)向驗(yàn)證曲線Fig. 5 Verification curves of calculated and measured wind direction at ocean stations

3.2 大區(qū)域潮流場(chǎng)模型潮位驗(yàn)證

本文搜集赤灣、珠海、大萬(wàn)山及廣州四個(gè)海洋站2017 年8 月（1713 號(hào)臺(tái)風(fēng)“天鴿”）及2018 年9月（1822 號(hào)臺(tái)風(fēng)“山竹”）非臺(tái)風(fēng)影響期間的潮位實(shí)測(cè)資料，與模型計(jì)算潮位進(jìn)行驗(yàn)證，結(jié)果表明，在各水文條件下，各潮位站模型與實(shí)測(cè)的潮位過(guò)程線吻合良好，模型的漲、落潮歷時(shí)和相位與實(shí)測(cè)資料一致（圖6）。

圖6 海洋站計(jì)算潮位和實(shí)測(cè)潮位驗(yàn)證曲線Fig. 6 Verification curves of calculated and measured tide level at ocean station

3.3 小區(qū)域模型潮汐、潮流、波浪驗(yàn)證

利用研究海域C1～C5 五個(gè)潮流測(cè)站和C1、C5 兩個(gè)臨時(shí)潮位站大潮期的實(shí)測(cè)潮流資料對(duì)模型計(jì)算得到的潮流進(jìn)行率定。C1、C5 潮位驗(yàn)證及C1～C5 測(cè)點(diǎn)流速流向驗(yàn)證結(jié)果分別見(jiàn)圖7 和圖8。

圖7 大潮期潮位驗(yàn)證曲線Fig. 7 Verification curves of the tide level during the spring tide

從潮位驗(yàn)證結(jié)果來(lái)看，大、小潮各潮位站計(jì)算值與實(shí)測(cè)值基本一致。計(jì)算的高、低潮位出現(xiàn)的時(shí)間與實(shí)測(cè)高、低潮位出現(xiàn)的時(shí)間都吻合得較好。潮位平均偏差在0.04 m 左右，誤差范圍小于10%，滿(mǎn)足JTS/T 231-2—2010《海岸與河口潮流泥沙模擬技術(shù)規(guī)程》（交通運(yùn)輸部水運(yùn)局, 2010）的要求。從流速、流向驗(yàn)證結(jié)果來(lái)看，大、小潮各測(cè)點(diǎn)流速、流向計(jì)算值與實(shí)測(cè)值吻合較好，相位偏差較小，流速過(guò)程與現(xiàn)場(chǎng)基本一致，平均流速計(jì)算值與實(shí)測(cè)值偏差在10%左右，平均流向小于10°，也滿(mǎn)足《海岸與河口潮流泥沙模擬技術(shù)規(guī)程》要求。因此，本模型計(jì)算得出的結(jié)果是合理的，模擬的流場(chǎng)反映了工程海域的潮流運(yùn)動(dòng)特征。

4 結(jié)果與分析

4.1 采砂區(qū)沖淤特征

自2018 年發(fā)現(xiàn)采砂坑以來(lái)，共開(kāi)展了歷時(shí)3 年（2019 至2021 年）的 8 次采砂坑及管線區(qū)水深地形測(cè)量和沉積物取樣工作?，F(xiàn)場(chǎng)調(diào)查資料充分，詳細(xì)查清了采砂區(qū)的坑內(nèi)淤積、邊界后退、邊坡變化情況，能夠準(zhǔn)確反映調(diào)查期間采砂坑沖淤發(fā)育趨勢(shì)和規(guī)律。8 次水深地形調(diào)查結(jié)果如圖9 所示。

為評(píng)估采砂坑的穩(wěn)定性和地形演化情況，本文分析了KP4 和KP5 兩個(gè)典型斷面（斷面位置見(jiàn)圖1）自2019 年3 月到2022 年9 月的坡面角度變化及坡頂高程變化，結(jié)果如表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

對(duì)比圖9 不同時(shí)期的水深地形可以看出：

1）2019 年3 月至2022 年9 月期間，采砂坑總體形態(tài)保持一致，未發(fā)生驟沖、驟淤現(xiàn)象。

2）采砂坑內(nèi)總體為淤積環(huán)境，呈現(xiàn)淤平趨勢(shì)，深坑淤積強(qiáng)度最大，其中北部深坑在2020 年4月前回淤明顯，2021 年至今無(wú)明顯沖淤現(xiàn)象；南部深坑總體為淤積現(xiàn)象，淤積變化不大。

3）2020 年9 月至2021 年1 月采砂坑北部深坑西側(cè)、采砂坑中部發(fā)生局部變深，為采砂活動(dòng)所致；2021 年1 月以后采砂坑內(nèi)未發(fā)現(xiàn)采砂活動(dòng)導(dǎo)致水深地形變化。

4）2020 年9 月前采砂區(qū)北部深坑淤積強(qiáng)度較大，與臨近時(shí)間采砂活動(dòng)關(guān)系密切。

整體看來(lái)，觀測(cè)期間采砂坑形態(tài)總體保持不變，坑內(nèi)淤積為主，無(wú)驟沖、驟淤現(xiàn)象，北部深坑2021 年1 月前強(qiáng)淤積過(guò)程與采砂活動(dòng)關(guān)系密切，2021 年1 月以后采砂坑整體淤積變化不大。

圖10 至圖12 分別表示采砂區(qū)及周邊海域1 a、2 a、5 a 后沖淤特征。通過(guò)對(duì)比分析，采砂坑的南北兩側(cè)及東側(cè)為沖刷環(huán)境，局部沖刷強(qiáng)度為0.5 m/a 左右；坑內(nèi)為淤積環(huán)境，局部最大淤積強(qiáng)度達(dá)到0.6 m/a 以上，整體淤積強(qiáng)度在0.2 m/a。長(zhǎng)期條件下，采砂區(qū)的南北兩側(cè)及東側(cè)為沖刷環(huán)境，局部沖刷深度分別約為1.1 m 和1.6 m，采砂區(qū)內(nèi)為淤積環(huán)境，整體淤積強(qiáng)度分別約為0.5 m/a 和1 m/a。采砂區(qū)北側(cè)邊界向西北方向擴(kuò)展，遠(yuǎn)離管道走向，對(duì)管道無(wú)影響。采砂區(qū)南側(cè)邊界向東南方向擴(kuò)展，與管道走向有一定重疊，導(dǎo)致銅鼓水道附近KP12～KP14 段管道沖刷，局部強(qiáng)度分別為0.3 m/a、0.8 m/a左右，未達(dá)到管道埋設(shè)，其他管道段仍呈現(xiàn)淤積環(huán)境。

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

圖11 現(xiàn)狀水深條件下2 a 后沖淤量分布Fig. 11 Distribution of scouring and silting amount under the case of present water depth after two years

圖12 現(xiàn)狀水深條件下5 a 后沖淤量分布Fig. 12 Distribution of scouring and silting amount under the case of present water depth after five years

4.2 采砂坑邊界變化

模擬結(jié)果（圖13）表明，采砂邊界變化主要發(fā)生在采砂坑南北邊界，東側(cè)靠近管線段采砂坑邊界變化略小。

圖13 模擬采砂坑邊界變化Fig. 13 Variations of the simulated boundary of the sand mining area

采砂坑南北邊界后退主要受控于漲落潮海流垂直沖刷。1 a、2 a 和5 a 后，采砂區(qū)北側(cè)邊界分別向西北方向呈環(huán)狀擴(kuò)展0.4、0.8 和1.2 km，采砂區(qū)南側(cè)邊界向東南方向呈環(huán)狀擴(kuò)展0.5、0.9 和1.5 km。采砂區(qū)東側(cè)邊界后退主要由凸體局部沖刷所致，1 a、2 a 和5 a 后向靠近管道方向擴(kuò)展6、11 和30 m，擴(kuò)展段為與管道KP4+120 至KP11+120之間平行的部分。

4.3 采砂坑邊坡變化

對(duì)于采砂區(qū)東側(cè)邊界的邊坡，選取KP8+387、KP10+387 處2 條典型斷面分析長(zhǎng)期沖刷條件下的采砂區(qū)邊坡變化特征，如圖14 所示。可以看出，采砂坑邊坡發(fā)育趨勢(shì)為坡底淤積、邊坡變緩、坡頂侵蝕。坡頂侵蝕會(huì)導(dǎo)致采砂坑邊界線向管線擴(kuò)展。

圖14 長(zhǎng)期沖刷作用下采砂區(qū)邊坡變化Fig. 14 Slope changes of the sand mining area under the action of long-term scouring

5 結(jié) 論

本文分析了深港海底管道路由區(qū)域2019 年3 月至2022 年9 月期間8 次采砂坑及管線區(qū)水深地形測(cè)量資料，通過(guò)數(shù)值模擬并與實(shí)測(cè)資料相互驗(yàn)證，掌握區(qū)域海區(qū)及采砂坑海域的動(dòng)力環(huán)境、沖淤規(guī)律特征。數(shù)值模擬方法能夠準(zhǔn)確反映調(diào)查期間采砂坑沖淤發(fā)育趨勢(shì)和規(guī)律，分析結(jié)果可為海底管道路由采砂區(qū)制定有效防護(hù)方案、采取有效防護(hù)措施提供科學(xué)依據(jù)，其主要結(jié)論如下。

1）采砂坑總體形態(tài)保持不變，未發(fā)生驟沖、驟淤現(xiàn)象。

2）采砂坑內(nèi)整體為淤積環(huán)境，采砂坑大坡度邊坡主要分布在東南部邊坡、南部邊坡，以及東北側(cè)局部；邊坡坡度總體維持不變。

3）采砂坑?xùn)|側(cè)邊緣走向與漲落潮主流向斜交，為沖刷環(huán)境，采砂坑邊緣突出部分進(jìn)行調(diào)整，形成凸體局部沖刷，采砂坑邊緣線整體后退。

4）模擬沖刷1 年后東側(cè)邊界整體后退約6 m，模擬沖刷2 a 后東側(cè)邊界整體后退約11 m，模擬沖刷5 a 后東側(cè)邊緣整體后退約30 m。