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

胡 慶 唐國樂 鄭川貴

（中電建路橋集團有限公司，北京 100048）

0 引言

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

1 工程概況

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

根據(jù)儲量核實報告，廣東省汕尾市管轄海域JH21-09 礦區(qū)海砂礦圈定回填用海砂礦體2 個，上部為“濱?！侵奚啊?、下部為“濱淺海砂”。采用地質塊段法估算，上部“濱?！侵奚啊钡V體控制資源量為578.24 萬m3，平均含泥量為16.03%；礦體推斷資源量488.55 萬m3，平均含泥量為15.57%；下部“濱淺海砂”礦體控制資源量為1387.69 萬m3，平均含泥量為16.07%；礦體推斷資源量為1052.99 萬m3，平均含泥量為16.26%。

2 三維水動力模型建立與驗證

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

2.1 三維水動力模型基本方程

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

2.1.1 連續(xù)方程

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

2.1.2 動量方程

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

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

2.1.3 初始條件

初始速度場、水位場均為0。

2.2 模型搭建及地形處理

2.2.1 模型搭建

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

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

圖1 三維潮流模型

2.3 模型驗證

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

表1 各個站點數(shù)據(jù)偏差/誤差

2.2.2 地形處理

進行數(shù)學模型計算時，根據(jù)項目區(qū)實測水下地形確定水深，以海平面為標準進行轉換計算，并將各個高程點賦值于模型節(jié)點中[4]。

2.2.3 邊界處理

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

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

3 模型計算結果與分析

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

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

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

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

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

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

3.2 施工后潮流場分析

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

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

圖3 施工后工程海域大潮漲急流場（中層）

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

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

4 結論

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