基于三維水動(dòng)力模型的海砂開采對(duì)海洋環(huán)境的影響分析

胡 慶 唐國樂 鄭川貴

（中電建路橋集團(tuán)有限公司，北京 100048）

0 引言

砂石資源在我國占據(jù)極其重要的地位，它不僅為基礎(chǔ)設(shè)施的建設(shè)提供了強(qiáng)大的支持，也為經(jīng)濟(jì)社會(huì)的發(fā)展提供了堅(jiān)實(shí)的保障。由于陸地砂資源逐漸減少，加上經(jīng)濟(jì)的快速增長，海洋砂資源變得越來越重要并且已成為當(dāng)前勘探開發(fā)的熱點(diǎn)[1]。但在海砂開采過程中，由于開采施工而影響海洋環(huán)境也是不可忽略的一大重點(diǎn)，因此采用技術(shù)手段模擬海砂開采后海洋環(huán)境變化并分析對(duì)海洋環(huán)境的影響，保證合理開采利用海砂，彌補(bǔ)陸地礦產(chǎn)資源的不足，發(fā)揮海砂資源的經(jīng)濟(jì)效益，保證國家重大項(xiàng)目的海砂供應(yīng)具有積極意義。

1 工程概況

海砂開采工程項(xiàng)目位于廣東省汕尾市管轄海域碣石灣施公寮島東南側(cè)，根據(jù)礦區(qū)鉆孔30 個(gè)組合樣碎屑礦物分析結(jié)果，海砂的礦物組成主要為石英，次要礦物為長石，其中石英含量為88.24%~95.54%。砂層呈松散結(jié)構(gòu)，平行層理構(gòu)造或塊狀構(gòu)造，緊密密度1700kg/m3，堆積密度1530kg/m3。

根據(jù)儲(chǔ)量核實(shí)報(bào)告，廣東省汕尾市管轄海域JH21-09 礦區(qū)海砂礦圈定回填用海砂礦體2 個(gè)，上部為“濱海—三角洲砂”、下部為“濱淺海砂”。采用地質(zhì)塊段法估算，上部“濱海—三角洲砂”礦體控制資源量為578.24 萬m3，平均含泥量為16.03%；礦體推斷資源量488.55 萬m3，平均含泥量為15.57%；下部“濱淺海砂”礦體控制資源量為1387.69 萬m3，平均含泥量為16.07%；礦體推斷資源量為1052.99 萬m3，平均含泥量為16.26%。

2 三維水動(dòng)力模型建立與驗(yàn)證

采砂區(qū)水深約15~20m，平均采砂深度為9.08m。工程所在海域以潮流作用為主，采砂施工使區(qū)域地形發(fā)生變化，潮流動(dòng)力受到一定影響，產(chǎn)生的懸浮泥沙隨海流輸運(yùn)擴(kuò)散，對(duì)采砂區(qū)及附近海域會(huì)產(chǎn)生一定影響。該文通過建立附近海域三維水動(dòng)力模型來分析采砂工程對(duì)海洋水文動(dòng)力環(huán)境的影響。

2.1 三維水動(dòng)力模型基本方程

三維水動(dòng)力數(shù)學(xué)模型是基于三維不可壓縮流，采用雷諾平均化N-S 方程，結(jié)合Boussinesq 假設(shè)和和靜水壓力假設(shè)，考慮湍流的影響、流體的密度變化，并考慮溫度平衡的因素，以更好地模擬水動(dòng)力學(xué)過程[2]。三維水動(dòng)力模型建立的基本方程主要包括連續(xù)方程（如公式（1）所示）和動(dòng)量方程（如公式（2）~公式（5）所示）。

2.1.1 連續(xù)方程

連續(xù)方程是流體力學(xué)的基本方程之一，反映的是流體運(yùn)動(dòng)和流體質(zhì)量分布的關(guān)系，是質(zhì)量守恒定律在流體力學(xué)中的應(yīng)用，具體如公式（1）所示。

2.1.2 動(dòng)量方程

動(dòng)量方程是流體運(yùn)動(dòng)的最基本運(yùn)動(dòng)學(xué)原理，即找出流體運(yùn)動(dòng)和它所受作用力之間的關(guān)系的數(shù)學(xué)表達(dá)式，依據(jù)的理論原理是牛頓的運(yùn)動(dòng)定律或動(dòng)量定理，具體如公式（2）~公式（5）所示。

式中：x、y、z為坐標(biāo)系3 個(gè)分量；D為總水深（m），D=H+η；H為平均海平面下的水深（m）；η為平均海平面起算水位（m）；u為x方向（東方向）流速（m/s）；v為y方向（北方向）流速（m/s）；w為σ坐標(biāo)系流速（m/s），；f為科氏參數(shù)；AM為垂直湍流黏滯系數(shù)；AH為水平湍流黏滯系數(shù)；ρ0為參考密度，取1025kg/m3；ρ為海水密度。

2.1.3 初始條件

初始速度場(chǎng)、水位場(chǎng)均為0。

2.2 模型搭建及地形處理

2.2.1 模型搭建

模型包括內(nèi)、外2 個(gè)模態(tài)。進(jìn)行計(jì)算時(shí)，外模態(tài)忽略垂直結(jié)構(gòu)，考慮水平對(duì)流和擴(kuò)散，計(jì)算二維變量；內(nèi)模態(tài)考慮垂向分層使用Sigma 坐標(biāo)，模型層數(shù)與實(shí)測(cè)海流對(duì)應(yīng)，分為3 層計(jì)算三維變量[3]。按穩(wěn)定性要求，水流數(shù)學(xué)模型的計(jì)算步長為30s。

模型求解采用非結(jié)構(gòu)網(wǎng)格中心網(wǎng)格有限體積求解，其優(yōu)點(diǎn)為計(jì)算速度較快。模型計(jì)算域涵蓋汕尾市、陸豐市海域，西至惠州雙月灣，東至陸豐石角尾山，海域東邊界為東經(jīng)114°50′19″，西邊界為東經(jīng)115°54′42″，南邊界至北緯22°19′43″，北至汕尾陸豐各個(gè)岸線。從外海至近海岸區(qū)域網(wǎng)格尺寸逐漸變小，其中近海岸海域網(wǎng)格尺寸為100~300m，外海網(wǎng)格尺度為2km~3km，采砂區(qū)海域周邊3km 范圍內(nèi)網(wǎng)格疏密根據(jù)需要確定。工程海域附近網(wǎng)格布置較密，非工程海域相對(duì)較疏，采砂區(qū)范圍內(nèi)網(wǎng)格尺寸為30m，采砂區(qū)周邊10km 范圍內(nèi)網(wǎng)格尺寸為30~100m。模型網(wǎng)格節(jié)點(diǎn)數(shù)為275276，單元數(shù)為111145，模型如圖1 所示。

圖1 三維潮流模型

2.3 模型驗(yàn)證

利用工程地區(qū)已有的4 個(gè)潮位站（C1~C4）的潮位資料和12 個(gè)流速測(cè)點(diǎn)（S1~S12）的流速、流向資料進(jìn)行模型的驗(yàn)證，并將模型計(jì)算所得的潮位、流速、流向數(shù)據(jù)與當(dāng)?shù)馗髡居涗浀膶?shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比，各個(gè)站點(diǎn)數(shù)據(jù)偏差/誤差見表1。

表1 各個(gè)站點(diǎn)數(shù)據(jù)偏差/誤差

2.2.2 地形處理

進(jìn)行數(shù)學(xué)模型計(jì)算時(shí)，根據(jù)項(xiàng)目區(qū)實(shí)測(cè)水下地形確定水深，以海平面為標(biāo)準(zhǔn)進(jìn)行轉(zhuǎn)換計(jì)算，并將各個(gè)高程點(diǎn)賦值于模型節(jié)點(diǎn)中[4]。

2.2.3 邊界處理

模型共設(shè)一個(gè)潮位開邊界，外海開邊界潮位以9 個(gè)調(diào)和常數(shù)的形式給出，由中國海域潮汐預(yù)報(bào)軟件China tide 計(jì)算獲得，主要考慮4 個(gè)半日分潮（M、S2、N 和Kz）、4 個(gè)全日分潮（Ki、Oi、Pr 和QI）及一個(gè)長周期分潮（Sa）[5]。

由表1 可以看出，模型計(jì)算得出的數(shù)據(jù)與當(dāng)?shù)赜涗浀膶?shí)測(cè)數(shù)據(jù)基本吻合，工程海域12 個(gè)潮流點(diǎn)的計(jì)算流速、流向和實(shí)測(cè)值較吻合，流速值的相對(duì)誤差都在8%以下，表明所建模型可以用來模擬研究工程實(shí)施造成的水動(dòng)力變化情況。

3 模型計(jì)算結(jié)果與分析

3.1 現(xiàn)狀潮流場(chǎng)分析

采用經(jīng)過驗(yàn)證的潮流數(shù)學(xué)模型，分別對(duì)該工程近水海域的大潮漲急、落急的表層、中層、底層潮流場(chǎng)進(jìn)行計(jì)算，結(jié)果如圖2 所示（以中層大潮漲急為例進(jìn)行展示）。

圖2 現(xiàn)狀工程海域大潮漲急流場(chǎng)（中層）

根據(jù)計(jì)算結(jié)果，漲潮時(shí)，海水從西邊界進(jìn)入碣石灣海域，并向西北向傳播。在碣石灣灣口和灣內(nèi)，受地形影響，潮流呈往復(fù)流。落潮時(shí)碣石灣內(nèi)海水朝東南方向流向外海。在施公寮島外圍形成回流，漲潮潮流呈逆時(shí)針方向流動(dòng)，落潮潮流呈順時(shí)針方向流動(dòng)。

采砂區(qū)位于施公寮島東南側(cè)約3.2km 處，大潮漲急流向大致在0°~45°，大潮落急流向大致在180°~225°。漲潮流流經(jīng)采砂區(qū)后傳播進(jìn)入施公寮島周邊，落潮流流經(jīng)該采砂區(qū)后流向外海。受島嶼地形影響，漲落潮在岬灣附近形成流速高值區(qū)。

工程區(qū)域水深較大，大潮漲潮平均流速在0.1~0.44m/s，大潮落潮平均流速在0.1~0.61m/s。工程海域潮流流速總體不大，大潮漲落潮平均流速在0.12m/s 左右，項(xiàng)目海區(qū)潮流流速整體較弱。

3.2 施工后潮流場(chǎng)分析

根據(jù)海域海砂資源儲(chǔ)量分析，按最不利情況考慮，即海砂開采區(qū)被平均挖深到24.08m。通過改變?cè)摬缮皡^(qū)海底地形，在采砂區(qū)周邊水域選取43 個(gè)代表點(diǎn)對(duì)施工后該工程近水海域的大潮漲急、落急的表層、中層、底層潮流場(chǎng)進(jìn)行計(jì)算，計(jì)算結(jié)果如圖3所示（以中層大潮漲急為例進(jìn)行展示），對(duì)比結(jié)果見表2（以海砂開采區(qū)中層大潮漲潮為例，選取JH21-09 區(qū)塊附近10 個(gè)點(diǎn)位進(jìn)行展示）。

表2 施工前、后中層大潮漲潮平均流速流向變化對(duì)比表

圖3 施工后工程海域大潮漲急流場(chǎng)（中層）

根據(jù)計(jì)算結(jié)果，采砂工程施工后，施工范圍海域的潮流影響相對(duì)較大并且表現(xiàn)為隨距離施工區(qū)域越遠(yuǎn)，影響越小的規(guī)律。采砂施工完成后，對(duì)施工海域南北方向的影響主要為潮流減速，對(duì)東西方向的影響主要為潮流加速，影響范圍為南北向和東西向約6km，工程前、后流速變化幅度均在0.022m/s 以下。

從工程區(qū)域大范圍海域流速變化來看，漲落潮時(shí)采砂工程對(duì)周邊海域產(chǎn)生的流速變幅除了采砂區(qū)內(nèi)及采砂區(qū)附近的t1、t18、t23、t27、t28、t31、t32、t36、t40 點(diǎn)位外，基本都在0.002m/s 及以下，流向變幅基本在1°以內(nèi)。t23、t28 位于工程區(qū)內(nèi)，流速變幅率在-15.9%~-21.4%，流向變化幅度在-1.1°~1.5°。t18、t32 靠近工程區(qū)，其流速變幅率不超過11%，流向改變幅度較小，不超過5°。因此可認(rèn)為采砂工程對(duì)采砂區(qū)內(nèi)流場(chǎng)的影響主要集中于采砂區(qū)內(nèi)，采砂區(qū)2.5km以外海域潮流場(chǎng)基本無變化，對(duì)大范圍海域潮流場(chǎng)的影響較小。

4 結(jié)論

該文通過采用數(shù)學(xué)模型對(duì)JH21-09 區(qū)塊所在海域的采砂工程進(jìn)行海洋環(huán)境影響模擬分析，得出如下結(jié)論：1）采砂工程施工完成后對(duì)海域水動(dòng)力條件的影響主要集中在海砂開采區(qū)2.5km 范圍內(nèi)，在2.5km 范圍外流速變化小于0.01m/s，流向變化小于1°。2）采砂工程對(duì)大范圍海域流場(chǎng)影響較小，對(duì)紅海灣開發(fā)區(qū)周邊海域不會(huì)造成明顯影響，不會(huì)造成海岸線坍塌等危害并且采砂工期短，結(jié)束后可恢復(fù)正常。