砂壩潟湖海灣水交換能力的數(shù)值模擬研究
——以茂名博賀灣為例

陳宏洺，張淑華，何 杰

(1.河海大學(xué) 海岸災(zāi)害及防護(hù)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210098；2.河海大學(xué) 港口海岸與近海工程學(xué)院， 南京 210098；3.南京水利科學(xué)研究院，南京 210029)

隨著經(jīng)濟(jì)發(fā)展和漁港漁業(yè)、工業(yè)碼頭迅速發(fā)展，導(dǎo)致博賀灣水體富營養(yǎng)化上升，未來港區(qū)規(guī)劃改變近海岸線布置。人類活動(dòng)、圍海造地、工業(yè)碼頭用地和海岸漁業(yè)都將影響海灣的有效水域面積和港灣水交換情況，對(duì)海灣生態(tài)系統(tǒng)、港口和航運(yùn)等有重大影響。水交換能力是評(píng)價(jià)海灣環(huán)境容量和環(huán)境質(zhì)量的重要指標(biāo)[1]，交換能力的強(qiáng)弱直接關(guān)系到海灣的水質(zhì)狀況。

目前已有多名國內(nèi)外學(xué)者開展了對(duì)河口、港灣水域的水交換研究，曾相明[2]等基于POM模型的基礎(chǔ)上,計(jì)算了不同時(shí)期的納潮量和水交換能力，得到較長時(shí)間段內(nèi)實(shí)施的大量圍填海工程納潮量和水交換累積影響明顯；陳振華[3]等在POM模型的基礎(chǔ)上建立了水質(zhì)模型，模擬欽州灣的水交換過程，結(jié)果表明欽州灣水交換能力整體上較強(qiáng)；張學(xué)慶[4]等基于EFDC模型,利用水齡概念,研究潮和風(fēng)對(duì)海州灣及毗鄰海域水交換的影響；姚姍姍[5]等通過模擬示蹤劑的濃度變化，模擬了海棠灣規(guī)劃在不同工況下的水體交換率，結(jié)果表明連續(xù)潮作用30 d后海棠灣水系整體水體交換率均在69%左右；Ren[6]等建立三維水齡模型，模擬了不同動(dòng)力條件下珠江口的水交換過程。何杰[7]等數(shù)學(xué)模型模擬了珠江口南沙港區(qū)挖入式港池的水體交換，模擬結(jié)果表明港池尾部水體交換率最低。董禮[8]先以溶解態(tài)的保守性物質(zhì)作為灣內(nèi)水的示蹤劑，建立對(duì)流-擴(kuò)散性的海灣水交換數(shù)學(xué)模型。戴志軍[9]等利用Gorden關(guān)系式，從茂名海域相衛(wèi)星影像定量反演了表層懸浮泥沙濃度的時(shí)空分布特征，結(jié)果表明茂名沿海水體表層懸浮泥沙濃度偏低,一般低于200 mg/L。

對(duì)博賀灣海域的潮流流場(chǎng)、水交換能力的科學(xué)認(rèn)識(shí)是對(duì)其科學(xué)合理開發(fā)的重要基礎(chǔ)，開展評(píng)估博賀灣污染物輸運(yùn)及水交換能力具有重要的科學(xué)意義。博賀灣為單口門海灣，較易存在水交換不暢的問題，國內(nèi)外針對(duì)博賀灣海域水交換能力的研究相對(duì)較少。本文建立了覆蓋博賀灣及其附近水域的平面二維潮流數(shù)學(xué)模型，在對(duì)近期實(shí)測(cè)水文資料驗(yàn)證的基礎(chǔ)上，模擬了博賀灣水交換過程及其動(dòng)力機(jī)制，分析了海灣的潮流、納潮量，文中采用一個(gè)大潮型作為分析博賀灣水體交換。

1 博賀灣自然條件

博賀灣位于廣東省西部茂名市電白區(qū)境內(nèi)，東與陽江市的海岸相接，西與水東港砂壩潟湖海灣毗連，當(dāng)?shù)厣钏x岸較近，潮汐動(dòng)力不強(qiáng)，陸源來沙極少，砂質(zhì)海床覆蓋大部分海區(qū)，水體含沙量很低，是廣東省優(yōu)質(zhì)海灣。博賀灣是一個(gè)典型的砂壩潟湖海灣，其地貌單元包括西側(cè)尖崗嶺—博賀的大砂壩和東側(cè)東閣嶺—蓮頭嶺的連島砂壩、面積約36 km2的潟湖、指向西南方向槽長約2 km的潮汐通道以及口門外側(cè)由砂咀和攔門沙淺灘組成的的水下落潮三角洲。博賀灣屬弱潮海區(qū)，平均潮差為1.05～2.80 m，最大潮差3.24 m，流運(yùn)動(dòng)形式總體上為往復(fù)流，受海岸地形影響，蓮頭嶺以東海域漲潮流向指向西偏南，落潮流向指向東偏北，繞過蓮頭嶺后漲潮流向指向西偏北，落潮流向指向東偏南。博賀灣內(nèi)存在大片淺灘，落潮期間會(huì)露出水面。

圖1為博賀灣水深及潮位、潮流觀測(cè)站位分布圖，T1,T2,T3為潮位測(cè)站，主要位于港池內(nèi)，H1～H6為潮流觀測(cè)站，分布為博賀灣外。博賀灣的潮流屬不正規(guī)半日潮流，為略帶旋轉(zhuǎn)的往復(fù)流運(yùn)動(dòng)，有明顯的駐波特征，漲落潮轉(zhuǎn)流出現(xiàn)在高、低潮平潮時(shí)刻，漲潮流歷時(shí)普遍大于落潮流歷時(shí)，漲潮平均流速一般均要大于落潮平均流速，2018年4月底水文實(shí)測(cè)各潮流站潮垂線平均流矢量圖如圖2所示。茂名博賀灣由于沒有大河入海，陸域來沙較少，泥沙主要以岸灘供沙為主。2004年10月實(shí)測(cè)平均含沙量為0.052～0.106 kg/m3，2018年5月實(shí)測(cè)平均含沙量為0.00～0.034 kg/m3，海域?qū)崪y(cè)懸浮泥沙含量很小。

圖1 博賀灣水深及潮位、潮流觀測(cè)站位分布Fig.1 Water depth and location of tide and tidal current observation in the Bohe bay圖2 附近海域各潮流站實(shí)測(cè)流速矢量Fig.2 Observed velocity vectors at observation stations nearby sea area

2 平面二維潮流數(shù)學(xué)模型

2.1 控制方程

在笛卡爾直角坐標(biāo)系下，根據(jù)靜壓和勢(shì)流假定，在水域面積較大的開闊水域，平面尺度遠(yuǎn)大于垂直尺度。將平面二維水流運(yùn)動(dòng)和物質(zhì)輸運(yùn)方程寫成如下的向量表示形式

(1)

式中：U=(d,du,dv,ds)Td為全水深d=h+ζ(h為水平面以下水深；ζ為潮位)；

2.2 數(shù)值求解

數(shù)值求解采用有限體積法對(duì)水沙方程進(jìn)行離散求解，實(shí)質(zhì)就是以單元為對(duì)象進(jìn)行水量、動(dòng)量和沙量的平衡，計(jì)算結(jié)果精度高，能處理含間斷或陡梯度的流動(dòng)。為了計(jì)算上的方便，統(tǒng)一采用三角形單元對(duì)計(jì)算區(qū)域進(jìn)行離散，并將單一的網(wǎng)格單元作為控制元，物理變量配置在每個(gè)單元的中心。

將第i號(hào)控制元記為Ωi，在Ωi上對(duì)向量式的基本方程組(1)進(jìn)行積分，并利用Green公式將面積分化為線積分，得

(2)

沿單元邊界線積分可以表示為三角形各邊積分之和

(3)

式(3)的求解主要分為三個(gè)部分，一是對(duì)流項(xiàng)的數(shù)值通量求解，二是紊動(dòng)項(xiàng)的求解，三是源項(xiàng)中底坡項(xiàng)的處理。對(duì)流項(xiàng)基面數(shù)值通量的求解格式有多種，這里采用Roe格式的近似Riemann解[10]。淺水方程的紊動(dòng)粘性項(xiàng)采用單元交界面的平均值進(jìn)行估算[11]，底坡源項(xiàng)采用特征分解法處理[12]。

2.3 水交換數(shù)值模擬

利用面源示蹤劑為保守性物質(zhì)，建立博賀灣水交換的數(shù)值模式，在計(jì)算過程中不考慮總投放物質(zhì)量衰減，隨水流運(yùn)動(dòng)、地形水深進(jìn)行擴(kuò)散和輸移。設(shè)海灣初始濃度場(chǎng)為C(x,y,d,t0)，瞬時(shí)濃度場(chǎng)為C(x,y,d,t)。則不同時(shí)刻，不同位置的灣內(nèi)水被外海水(示蹤物濃度為0)置換的比率R(x,y,d,t)為

(4)

若設(shè)C(x,y,d,t0)=1.0 mg/L，則(4)簡化為

R(x,y,d,t)=1-C(x,y,d,t)

(5)

3 網(wǎng)格剖分及模型驗(yàn)證

3.1 網(wǎng)格剖分及模型參數(shù)

圖3 網(wǎng)格剖分效果示意Fig.3 Splitted grids of model

茂名博賀新港西區(qū)處在博賀灣內(nèi)，東區(qū)處在蓮頭嶺以東。數(shù)學(xué)模型東邊界設(shè)在水東灣以西，西邊界位于沙扒港以東，南邊界設(shè)在蓮頭嶺以南20 km的外海一線。數(shù)學(xué)模型東西長度49 km，南北向30 km，控制范圍面積接近1 300 km2。計(jì)算區(qū)域共劃9萬余個(gè)三角形網(wǎng)格，對(duì)博賀灣外海域水面寬闊，采用較為稀疏的大尺度網(wǎng)格剖分，最大尺度為1 000 m(剖分效果見圖3)。對(duì)港區(qū)內(nèi)主要水域進(jìn)行網(wǎng)格加密，網(wǎng)格最小尺度為20 m。這種非結(jié)構(gòu)網(wǎng)格劃分法可以克服矩形網(wǎng)格鋸齒形邊界所造成的流動(dòng)失真，也可以避免生成有結(jié)構(gòu)貼體曲線網(wǎng)格的復(fù)雜計(jì)算和其他困難。模型計(jì)算參數(shù)如表1所列。

圖4 博賀灣水域示蹤劑布置范圍示意Fig.4 Schematic diagram of tracer layout range in Bohe bay

表1 模型計(jì)算參數(shù)Tab.1 Model calculation parameters

3.2 示蹤劑布置

模擬在博賀灣和博賀漁港布置面源示蹤劑，水域面積34.2 km2，初始濃度為1.0，在潮流作用下，與外海示蹤劑為0的干凈水交換，布置示意圖如圖4。模型計(jì)算結(jié)果的濃度值為示蹤劑指標(biāo)，即原始物質(zhì)經(jīng)過擴(kuò)散輸移過程后的殘留度。博賀灣內(nèi)水深較淺，落潮期間灣內(nèi)有大片淺灘露出水面，面源示蹤劑將隨水體流出博賀灣，落潮期間會(huì)露出水面，由于數(shù)模中設(shè)置了最小閾值水深來控制干濕，這種方式將導(dǎo)致物質(zhì)濃度的不守恒，誤差對(duì)于結(jié)果影響較小。即假定示蹤劑完全隨水流運(yùn)動(dòng)，在示蹤劑流出博賀灣時(shí)在淺灘海床上不會(huì)有殘余粘附。

3.3 模型驗(yàn)證

潮流泥沙數(shù)學(xué)模型與天然相似的條件，主要取決于模型計(jì)算出的潮流場(chǎng)和含沙量場(chǎng)與實(shí)測(cè)結(jié)果的吻合程度。通過2018年4月底的夏季大潮水文泥沙實(shí)測(cè)資料對(duì)數(shù)學(xué)模型進(jìn)行率定。該次水情組合有3站潮位和6條垂線的流速、流向資料可供模型驗(yàn)證，潮位站位置分布見圖1。

圖5為模型采用夏季大潮水情對(duì)潮位站實(shí)測(cè)潮位、垂線流速和流向過程的驗(yàn)證。圖5中列舉了T1,T2,H1,H2潮位驗(yàn)證過程，模型得到潮位過程與天然情況基本吻合，模擬各條垂線的流速、流向過程模擬與實(shí)測(cè)趨于相似。計(jì)算得到個(gè)別站位外其他站的高、低潮位計(jì)算值與實(shí)測(cè)值誤差均在±0.10 m范圍以內(nèi)；垂線漲潮、落潮平均流速的計(jì)算值與實(shí)測(cè)值相差均在10%以內(nèi)；驗(yàn)證結(jié)果總體符合數(shù)學(xué)模型的精度要求。

4 博賀灣水動(dòng)力分析

4.1 流態(tài)變化

圖6中列舉了現(xiàn)狀條件下水域平面漲落潮流場(chǎng)分布，漲落潮的平均流速是垂向平均流速在整個(gè)漲潮時(shí)間內(nèi)或落潮時(shí)間內(nèi)的流速平均值。由圖6可知在灣口水域形成漲潮流由西北轉(zhuǎn)向東北、落潮流由西北轉(zhuǎn)向西南的主流流態(tài)，其因?yàn)樯忣^半島類似一條巨大的鉤頭丁壩，使進(jìn)出博賀灣的潮流繞著壩頭作貼岸運(yùn)動(dòng)。博賀灣入口處深槽水域流速較大，漲落潮最大流速可分別達(dá)到0.8 m/s左右；博賀灣內(nèi)淺灘水域流速相對(duì)較弱，流速普遍在0.2 m/s左右；落潮過程中，落潮流沿灣內(nèi)潮溝向?yàn)惩獬隽?，博賀灣內(nèi)有大片淺灘露出水面。

5-a T1潮位過程驗(yàn)證5-b T2潮位過程驗(yàn)證5-c H1流速過程驗(yàn)證

5-d H2流速過程驗(yàn)證5-e H1流向過程驗(yàn)證5-f H2流向過程驗(yàn)證圖5 潮位、流速和流向過程驗(yàn)證Fig 5 Verification of tide level, flow rate and flow direction

6-a 漲急時(shí)刻 6-b 落急時(shí)刻圖6 博賀灣水域平面流態(tài)Fig.6 Flow pattern of Bohe bay waters

4.2 水交換能力

水交換能力反映了一個(gè)半封閉海灣接受外來“新水”的能力，是衡量海灣自凈能力大小的重要指標(biāo)。2018年5月份夏季大潮水情時(shí)為例，圖7顯示的是現(xiàn)狀條件下蹤劑隨時(shí)間分布典型情況可知，投放初期在一個(gè)潮周期的高潮位；在投放后6 h達(dá)到低潮位，落潮期間有大片淺灘露出水面，絕大部分示蹤劑隨深槽集中在灣外，在蓮頭嶺北側(cè)灣澳處有殘留；在投放后13 h經(jīng)歷下一個(gè)漲潮，水體大部分隨著漲潮流重新進(jìn)入博賀灣，部分水體被蓮頭嶺隔開后向南流；在投放后19 h，再次退潮水體絕大部分離開博賀灣，灣外示蹤劑濃度相比6 h減弱，部分水體繞過蓮頭嶺向東南方向拓展。博賀灣水體示蹤劑指標(biāo)值隨漲、落潮變化，落潮期間指標(biāo)值下降，漲潮期間指標(biāo)值上升，博賀灣水體示蹤劑指標(biāo)值整體上呈減小趨勢(shì)。示蹤劑在漲、落潮過程中不斷進(jìn)入和退出灣口，博賀灣內(nèi)示蹤劑的指標(biāo)在不斷下降，部分示蹤劑隨東南向落潮流進(jìn)入到蓮頭嶺以南水域。表明博賀灣具有良好的水體交換能力。

7-a 投放初期7-b 投放后6 h

7-c 投放后13 h7-d 投放后19 h圖7 博賀灣示蹤劑隨時(shí)間分布情況Fig.7 Tracer distribution over time in Bohe bay

8-a 投放初期8-b 投放后6 h

8-c 投放后13 h8-d 投放后19 h圖8 粵西沿岸流作用示蹤劑分布情況Fig.8 Tracer distribution along the coastal flow in the West Coast of Guangdong Province

4.3 沿岸流對(duì)水交換影響

粵西近岸地區(qū)常年均存在一股向西流動(dòng)的沿岸海流，稱為“粵西沿岸流”。粵西沿岸流由于受到風(fēng)及外海海流格局的影響，有顯著的季節(jié)變化；楊毅等[13]分析了粵西沿岸海流特征，粵西沿岸海流主要為西南向流，只有在夏季西南季風(fēng)較大時(shí)才會(huì)有東北向流；楊士瑛等[14]研究發(fā)現(xiàn)粵西沿岸流對(duì)徑流、污染物等運(yùn)移和擴(kuò)散起著決定性的作用；本次模擬了夏季一個(gè)潮周期在沿岸流作用下的博賀灣典型時(shí)刻的水交換變化情況。

如圖8所示，分別選取了4個(gè)典型時(shí)刻可以看出漲潮流繞過蓮頭嶺東北向進(jìn)入博賀灣，落潮流則西南向流出博賀灣后，沿岸線方向繼續(xù)向西南流動(dòng)。由于落潮流方向的改變，沿岸流將會(huì)一定程度上對(duì)博賀灣水交換影響。博賀灣內(nèi)示蹤劑隨落潮流流出灣口，在沿岸流作用下，示蹤劑偏向西側(cè)流動(dòng)，對(duì)比圖7可知在沿岸流作用下水體交換能力有一定的增強(qiáng)；漲潮期間，部分示蹤劑隨張潮流進(jìn)入博賀灣，仍有一部分殘留在水東灣與博賀灣之間的近岸水域；再次退潮過程中，大部分示蹤劑已完全流出博賀灣，且向西移動(dòng)范圍超過了水東灣，說明了粵西沿岸海一定程度上加強(qiáng)了博賀灣水體交換的能力。

5 博賀灣規(guī)劃對(duì)水交換影響

根據(jù)茂名博賀新港區(qū)總體規(guī)劃[15]，未來博賀港東側(cè)白沙頭至南側(cè)蓮頭嶺沙咀沿岸均規(guī)劃有填海項(xiàng)目。博賀灣內(nèi)，采用順岸布置，主要規(guī)劃裝備制造等臨港產(chǎn)業(yè)，博賀灣規(guī)劃完成后，博賀灣內(nèi)填海面積為7.8 km2，航道、港池開挖面積為3.9 km2，開挖底標(biāo)高為-10.0 m，進(jìn)出港航道的航槽底寬為200 m。小面積填海工程對(duì)海灣納潮量和流場(chǎng)的影響不大,但較長時(shí)間段內(nèi)多個(gè)圍填海工程累積起來的影響卻比較明顯，本文探討了博賀灣港東、西區(qū)的各類碼頭規(guī)劃完工后對(duì)灣內(nèi)水動(dòng)力的影響。

表2 博賀灣各斷面納潮量統(tǒng)計(jì)Tab.2 Tidal prism in each section of Bohe bay

5.1 納潮量變化

由于海洋的潮汐現(xiàn)象，一個(gè)潮周期內(nèi)把從低潮到高潮海灣所能容納海水的量稱為納潮量，對(duì)海灣內(nèi)的生態(tài)、環(huán)境及沖淤變化等方面的意義重大。納潮量的大小可以直接影響海灣與外海的交換程度,從而制約海灣的自凈能力，因此對(duì)維持海灣的良好生態(tài)環(huán)境至關(guān)重要。由表2統(tǒng)計(jì)了博賀灣和博賀漁港兩個(gè)斷面(位置見圖4所示)在現(xiàn)狀、規(guī)劃方案實(shí)施后的納潮量及其變化幅度。通過對(duì)比可以看出：在低潮位博賀灣中水體潮量為14.9億m3，在高潮位潮量為64.1億m3，即一個(gè)漲潮過程中，博賀灣有78.78%的新水涌入博賀灣中；在規(guī)劃實(shí)施后，DM1斷面納潮量減小0.42%，博賀漁港及其內(nèi)灣納潮量基本不受博賀新港規(guī)劃工程的影響；博賀灣斷面DM2潮量將減小16.4%，這是由于博賀新港規(guī)劃工程西區(qū)陸域圍填面積達(dá)到7.8 km2，圍填范圍超出蓮頭嶺的西北角，影響水體進(jìn)出港灣。

5.2 水交換變化

模擬灣內(nèi)初始示蹤劑濃度為1.0，通過數(shù)學(xué)模型模擬半月連續(xù)潮的水體交換過程，模型計(jì)算過程中的潮流控制邊界條件，如圖9所示采用15 d實(shí)測(cè)潮位資料進(jìn)行分析，連續(xù)半月潮中包括若干大潮、中潮、小潮，模擬時(shí)間為15 d，分別模擬了兩種工況在博賀灣水體交換情況。

由圖10可對(duì)比兩種工況在博賀灣水體示蹤劑指標(biāo)，可以看出：由于博賀灣規(guī)劃實(shí)施后水域面積有所減小，港池水域示蹤劑在落潮期間仍留在原處，漲潮期間示蹤劑隨漲潮流重新進(jìn)入博賀灣，落潮期間港池水域的示蹤劑比較難以隨落潮流流出灣口，因此，港池水域的水體交換能力相對(duì)要差一些；水交換速率在示蹤劑投放初期會(huì)比現(xiàn)狀條件弱一些，但后期水交換情況基本與現(xiàn)狀一致；現(xiàn)狀條件下一個(gè)全日潮大潮水情的水體交換率達(dá)到69%，中潮水情下水體交換率達(dá)到55%，小潮水情下水體交換率達(dá)到43%，由于大小潮的不同潮位過程、不同潮汐動(dòng)力造成水體交換率隨著潮汐變化而變化；博賀灣達(dá)到90%水交換的時(shí)間為3～4 d，達(dá)到50%水交換時(shí)間為0.5 d，進(jìn)一步說明了博賀灣具有良好的水體交換能力。

圖9 博賀灣潮位過程線Fig.9 Tide level process line in Bohe bay圖10 博賀灣水交換時(shí)間的示蹤劑指標(biāo)變化Fig.10 Change of total tracer during the water exchange process in Bohe bay

6 結(jié)論

本文建立了博賀灣二維水質(zhì)模型，通過與實(shí)測(cè)資料對(duì)比，證明了模型的可靠性。利用模型研究分析了規(guī)劃對(duì)于博賀灣水動(dòng)力環(huán)境的影響，可得：

(1)現(xiàn)狀條件下博賀灣內(nèi)存在大片淺灘，落潮期間會(huì)露出水面，灣內(nèi)布置的面源示蹤劑大部分會(huì)隨落潮流到灣外，大潮水情下一個(gè)全日潮過程灣內(nèi)水體交換率達(dá)到69%，中潮水情水體交換率達(dá)到55%，小潮水情水體交換率達(dá)到43%，博賀灣總體上水交換情況良好。

(2)博賀新港規(guī)劃階段實(shí)施后，陸域圍填面積達(dá)到7.8 km2，灣口斷面納潮量將減小16.4%；博賀漁港及其內(nèi)灣納潮量基本不受博賀新港規(guī)劃工程的影響。

(3)規(guī)劃實(shí)施后，水交換速率在示蹤劑投放初期會(huì)比現(xiàn)狀條件弱一些，后期水交換情況基本與現(xiàn)狀一致，這也進(jìn)一步說明了博賀灣具有良好的水體交換能力。