珠江出?？谙坛睌?shù)值模擬技術(shù)

侯堋，舒勰俊

(1.珠江水利科學(xué)研究院 廣州 510611；2.國(guó)家海洋局海洋咨詢中心 北京 100860；3.國(guó)家海洋局南海規(guī)劃與環(huán)境研究院 廣州 510300)

?

珠江出?？谙坛睌?shù)值模擬技術(shù)

侯堋1，舒勰俊2，3

(1.珠江水利科學(xué)研究院 廣州 510611；2.國(guó)家海洋局海洋咨詢中心 北京 100860；3.國(guó)家海洋局南海規(guī)劃與環(huán)境研究院 廣州 510300)

珠江三角洲是中國(guó)經(jīng)濟(jì)發(fā)展速度最快的地區(qū)之一，但珠江出海口咸潮上溯的情況日益嚴(yán)重，已經(jīng)嚴(yán)重影響珠江口的取水安全；對(duì)每年枯水期珠江口咸潮上溯的情況進(jìn)行預(yù)警預(yù)報(bào)，從而采取相應(yīng)的抑咸措施具有重大意義。文章通過建立珠江三角洲的二、三維嵌套水流鹽度數(shù)學(xué)模型，采用2005年1月的枯水期實(shí)測(cè)潮流鹽度資料對(duì)模型進(jìn)行率定，數(shù)學(xué)模型的計(jì)算結(jié)果與實(shí)測(cè)擬合較好，說明該模型運(yùn)用于珠三角咸潮數(shù)值模擬的可靠性，為咸潮預(yù)警預(yù)報(bào)技術(shù)的下一步研究奠定基礎(chǔ)。

珠江三角洲；咸潮入侵；海洋災(zāi)害；防災(zāi)減災(zāi)

珠江三角洲是中國(guó)經(jīng)濟(jì)發(fā)展速度最快的地區(qū)之一，水資源豐富，水系復(fù)雜，形成獨(dú)有的“三江匯流，八口入海”的流域特色。近幾年珠江三角洲干旱比較嚴(yán)重，北江、西江和東江的水流量普遍減少，尤其是西江的水流量減至自1903年以來的最低水平。2005年春季珠江三角洲暴發(fā)特大咸潮災(zāi)害，咸潮入侵導(dǎo)致江水鹽度升高，影響江水中營(yíng)養(yǎng)鹽的濃度與分布，造成水體生態(tài)環(huán)境改變，嚴(yán)重危害口門處居民的取水安全[1]。因此，采用數(shù)值模擬技術(shù)對(duì)珠江三角洲枯季咸潮上溯情況進(jìn)行預(yù)警預(yù)報(bào)，結(jié)合上游水情采取抑咸措施，確保生態(tài)環(huán)境質(zhì)量和口門取水安全，對(duì)促進(jìn)珠江三角洲地區(qū)及港澳特區(qū)的社會(huì)經(jīng)濟(jì)穩(wěn)定和持續(xù)發(fā)展都有重要意義[2]。

1 珠江三角洲二、三維水流鹽度模型的建立

對(duì)珠江三角洲進(jìn)行咸潮預(yù)警預(yù)報(bào)，需要建立一個(gè)范圍較大、能涵蓋各主要因素的水流鹽度數(shù)學(xué)模型，對(duì)珠江三角洲鹽水入侵進(jìn)行全方位的系統(tǒng)研究。上邊界須取在不受潮汐控制的恒定流區(qū)域，下邊界須取在外海-30 m等深線處。

二維模型模擬的區(qū)域?yàn)榫W(wǎng)河區(qū)，模型的上邊界為西江梧州，北江飛來峽、老鴉崗，東江三角洲大盛、麻涌、漳澎、泗盛圍，潭江石咀等水文站并采用實(shí)測(cè)流量過程，下邊界取至虎門大虎站、蕉門南沙站、洪奇門馮馬廟站、橫門橫門站、磨刀門燈籠山站、雞啼門黃金站、虎跳門西炮臺(tái)站和崖門官?zèng)_站。二維模型的網(wǎng)格數(shù)為154 992個(gè)。

三維模型模擬的區(qū)域?yàn)榘舜罂陂T和外海水域，上邊界為三角洲八大口門出口控制水文站(即網(wǎng)河區(qū)二維模型的下邊界)，下邊界取至外海-30 m等深線，采用荷包島、大萬山、擔(dān)桿頭等臺(tái)站的實(shí)測(cè)潮位資料。三維模型的平面網(wǎng)格數(shù)為33 111個(gè)，垂向分10層。二、三維模型的平面總網(wǎng)格數(shù)為188 103個(gè)。

1.1 水動(dòng)力數(shù)學(xué)模型

1.1.1 二維水動(dòng)力數(shù)學(xué)模型

深度平均的連續(xù)方程如下：

式中：Q為因水的流出和流進(jìn)、降雨和蒸發(fā)引起的每單位面積的流量，且

qin和qout分別是當(dāng)?shù)孛矿w積單元水源和匯(L/s)；P為降雨的非當(dāng)?shù)卦错?xiàng)；E為蒸發(fā)的非當(dāng)?shù)卦错?xiàng)。標(biāo)注這些吸入量，如電廠吸收水將被作為一個(gè)匯，在自由表面降雨是一個(gè)源項(xiàng)、蒸發(fā)是一個(gè)匯。式中U、V為、η方向流速分量，為水位分量，G為拉梅系數(shù)。

水平方向動(dòng)量方程

在ξ-和η-方向的動(dòng)量方程如下：

密度變化被忽略，除在斜壓項(xiàng)，Pξ和Pη代表壓力梯度；動(dòng)量方程中的Fξ和Fη代表水平雷諾應(yīng)力的不平衡；Mξ和Mη代表源于外部動(dòng)量的源和匯的貢獻(xiàn)(外部驅(qū)動(dòng)、水工建筑物、排放或水的提取、波浪應(yīng)力等)。

1.1.2 三維水動(dòng)力數(shù)學(xué)模型

三維水動(dòng)力數(shù)學(xué)模型必須在平面二維水動(dòng)力數(shù)學(xué)模型的基礎(chǔ)上考慮垂向變化。

σ-坐標(biāo)系下的垂直流速ω從連續(xù)方程中計(jì)算得：

在表層，降雨和蒸發(fā)的影響都考慮在內(nèi)。垂直流速ω定義在等σ-表面。ω是相對(duì)運(yùn)動(dòng)的σ-平面的垂直速度。ω是相對(duì)移動(dòng)的σ-平面的垂直流速，它解釋為同上升和沉降運(yùn)動(dòng)聯(lián)合的流速。在笛卡兒坐標(biāo)系的“物理”垂直流速w不包含在方程中，物理垂直流速的計(jì)算只在后處理用途中要求。這些流速可用水平流速、水深、水位和垂直ω-流速表達(dá)，依照：

1.2 鹽度數(shù)學(xué)模型

在模型中，鹽度輸運(yùn)使用一個(gè)三坐標(biāo)方向的對(duì)流-擴(kuò)散方程。源項(xiàng)和匯項(xiàng)包括模擬排放和收回。考慮一階衰變過程，一階衰變過程符合一個(gè)數(shù)值解法，指數(shù)衰減更復(fù)雜的過程。

輸運(yùn)方程在水平方向曲線網(wǎng)格坐標(biāo)和垂直方向σ坐標(biāo)，用一種守恒形式的公式表示：

λd(d+)c+S

λd代表第一階衰變過程和S是由水排放qin或收回qout，通過自由表面的熱交換Qtot引起的每單位面積的源項(xiàng)和匯項(xiàng)：

垂直湍流黏性系數(shù)DV定義為：

湍流閉合模型為解決所有其他無法解決的混合形式，便利地指定一個(gè)背景或者“周圍的”垂直混合系數(shù)。如，在強(qiáng)層結(jié)流，交界面的湍流渦擴(kuò)散率減少到0，垂直混合減少到分子擴(kuò)散。這個(gè)在本質(zhì)上不真實(shí)，因?yàn)閷a(chǎn)生內(nèi)部波。需指定一個(gè)背景混合參數(shù)Dback，如下：

二維深度平均模擬，水平渦擴(kuò)散DH也包含水平流垂直變化的貢獻(xiàn)(泰勒剪切彌散)，這部分不被深度平均3D對(duì)流擴(kuò)散方程考慮。深度平均輸運(yùn)模式的應(yīng)用，水體需很好的混合，泰勒剪切彌散次要，否則要使用三維模型。

1.3 二、三維模型嵌套技術(shù)

本文中珠江三角洲網(wǎng)河區(qū)采用二維平面模型、八大口門及外海處采用三維數(shù)學(xué)模型，因此需要采用嵌套技術(shù)將二維模型和三維模型進(jìn)行耦合。本文采用區(qū)域分裂法進(jìn)行聯(lián)解，區(qū)域分裂法是把待求的整個(gè)區(qū)域按不同的計(jì)算方法分成幾個(gè)子區(qū)域，在區(qū)域與區(qū)域交界之間選取公共節(jié)點(diǎn)，各子區(qū)域按給定的計(jì)算方法和計(jì)算步長(zhǎng)分別計(jì)算，計(jì)算時(shí)幾個(gè)子區(qū)域交替進(jìn)行，通過公共節(jié)點(diǎn)位相臨子區(qū)域提供邊界上的水力要素，以此將子區(qū)域聯(lián)接起來構(gòu)成一個(gè)整體的數(shù)學(xué)模型[3]。本文將計(jì)算區(qū)域分成兩個(gè)子區(qū)域，即二維模型計(jì)算部分和三維模型計(jì)算部分，三維模型的上邊界由二維模型的計(jì)算結(jié)果提供，二維子區(qū)域的下邊界條件由三維模型計(jì)算結(jié)果提供，子區(qū)域間計(jì)算因子相互傳遞，結(jié)合成一個(gè)整體。

2 三角洲二、三維水流鹽度模型驗(yàn)證

2.1 模型驗(yàn)證條件

(1)地形資料。模型驗(yàn)證地形資料采用1∶500 0的河道地形圖，其中西江干流(梧州至磨刀門)主要采用2005年在該地區(qū)測(cè)量的河道地形資料，其余八大口門采用最新測(cè)量的地形資料，其他地形采用1999年的實(shí)測(cè)地形資料。

(2)水文資料。模型驗(yàn)證水文資料取自近年來資料較為完整的史冊(cè)水文系列，選用2005年1月18日—2月4日的枯水期水文組合對(duì)模型進(jìn)行驗(yàn)證 。

2.2 模型驗(yàn)證成果

模型的計(jì)算時(shí)間為2005年1月18日9：00—2月4日12：00，通過實(shí)測(cè)資料對(duì)網(wǎng)河區(qū)及八大口門的采樣點(diǎn)進(jìn)行驗(yàn)證，由于篇幅關(guān)系，文中僅列出八大口門觀測(cè)點(diǎn)的潮位驗(yàn)證圖和4個(gè)采樣點(diǎn)流速驗(yàn)證圖，同時(shí)還列舉3個(gè)采樣點(diǎn)的鹽度驗(yàn)證圖(圖1～圖3)。

圖1 八大口門觀測(cè)點(diǎn)潮位驗(yàn)證

圖2 觀測(cè)點(diǎn)表底層流速驗(yàn)證

圖3 觀測(cè)點(diǎn)表底層鹽度驗(yàn)證

從驗(yàn)證成果可見：枯水水文條件下各潮位站模型與原型的潮位過程線吻合較好，模型的漲、落潮歷時(shí)和相位與原型實(shí)測(cè)資料基本一致，滿足精度要求；該水文條件下流速驗(yàn)證過程與實(shí)測(cè)過程吻合較好，相位基本一致，并且表層流速明顯大于底層流速，符合自然運(yùn)動(dòng)規(guī)律；鹽度的驗(yàn)證體現(xiàn)其隨潮流漲落而升降的特征，在網(wǎng)河區(qū)濃度由口門往上游沿程遞減，基本反映出鹽度在徑潮動(dòng)力共同作用下于西北江三角洲網(wǎng)河區(qū)的分布規(guī)律，即從垂向上看底層鹽度大于表層鹽度，并且磨刀門內(nèi)的站位在小潮時(shí)能夠形成第二次的鹽度小高峰，這與實(shí)測(cè)資料相符合。從垂向的流速和鹽度分布來看，垂向分布的情況基本符合現(xiàn)場(chǎng)資料分析結(jié)果，證明該模型能很好地復(fù)演珠江口鹽水入侵的三維特征。

3 結(jié)論

(1)構(gòu)建涵蓋珠江三角洲網(wǎng)河及河口區(qū)的二、三維耦合的水流鹽度數(shù)學(xué)模型，該模型在數(shù)值模擬方面解決因珠江三角洲網(wǎng)河數(shù)量多、交錯(cuò)復(fù)雜的水系模擬導(dǎo)致的計(jì)算耗時(shí)過長(zhǎng)的問題，同時(shí)在河口區(qū)采用三維數(shù)值模擬能更精確擬合咸水上溯過程中的垂向分布形式。

(2)通過選用2005年1月18日-2月4日的枯水期水文組合對(duì)模型進(jìn)行驗(yàn)證，結(jié)果表明本文的數(shù)學(xué)模型計(jì)算結(jié)果與實(shí)測(cè)資料吻合較好，能夠體現(xiàn)鹽水入侵的三維特征，證明數(shù)學(xué)模型的可靠性，為咸潮預(yù)警預(yù)報(bào)技術(shù)的下一步研究奠定基礎(chǔ)。

[1] 殷建平，王友紹.特大咸潮對(duì)珠江入海河段環(huán)境要素的影響[J].熱帶海洋學(xué)報(bào)，2006(4)：79-84.

[2] 徐峰俊，朱士康，劉俊勇.珠江河口區(qū)水環(huán)境整體數(shù)學(xué)模型研究[J].人民珠江，2003(5)：12-18.

[3] 謝作濤，羅景.長(zhǎng)江口一、二維嵌套水流鹽度數(shù)學(xué)模型[J].武漢大學(xué)學(xué)報(bào)：工學(xué)版，2007，40(2)：7-12.

The Numerical Simulation Model of Predicting the Salt Intrusion in the Pearl River Delta

HOU Peng1，SHU Xiejun2，3

(1.Pearl River Hydraulic Research Institute，Guangzhou 510611，China；2.Consultation Center of State Oceanic Administration，Beijing 100860，China；3.South China Sea Institute of Planning and Environmental Research，SOA，Guangzhou 510300，China)

The Pearl River Delta is one of the most developed area in china，but the situation of the saltwater intrusion has become more and more serious.The saltwater intrusion has influenced the safety of getting freshwater.It is necessary to predict the situation of the saltwater intrusion in dry seasons and take some action to restrain the intrusion.The paper introduced 2-D and 3-D coupling flow and salinity mathematical models to study the tidal current fields and salt data in January，2005.The relatively well fit result between computation output and the field data indicated that the model could reasonably simulate the hydrodynamic and salt fields，which could be the base of the research on the saltwater intrusion prediction technology.

The Pearl River Delta，Saltwater intrusion，Marine disaster，Disaster prevention and mitigation

2016-02-25；

2016-07-04

侯堋，碩士，研究方向?yàn)楹涌谂c海岸動(dòng)力學(xué)，電子信箱：houpen@126.com

P7

A

1005-9857(2016)08-0065-05