海平面上升對珠江河口東四口門鹽水入侵的影響

劉忠輝，宏波

(華南理工大學(xué)，廣東廣州510641)

鹽水入侵是由于河口鄰近海域的高鹽水團(tuán)隨潮汐漲潮流沿著潮汐通道向上游推進(jìn)，咸淡水混合使河口上游區(qū)域水體變咸的現(xiàn)象，會對沿途區(qū)域的水環(huán)境、農(nóng)業(yè)灌溉和生活供水等產(chǎn)生重要影響。珠江三角洲城市群擁有城鎮(zhèn)人口5 500萬，已成為全世界最大的城市片區(qū)[1]。該區(qū)域約占中國國土面積0.5%的土地，但其生產(chǎn)總值卻占中國總量的近20%[2]，是中國最重要的經(jīng)濟(jì)區(qū)域之一。在全球海平面上升的背景下，珠江三角洲地區(qū)已成為全球受鹽水入侵嚴(yán)重影響的地區(qū)之一[3]。從2000年開始，珠江三角洲沿海城市幾乎每年都會受到由于鹽水入侵導(dǎo)致的淡水供應(yīng)不足問題[4]。2005年12月，珠海、中山和廣州等地均受鹽水入侵影響，取水口水閘氯度嚴(yán)重超標(biāo)，最高達(dá)7 500 mg/L[5]。

眾多研究表明，海平面上升會強(qiáng)化鹽水入侵。Rice等[6]利用HEM-3D模型對美國東海岸切薩皮克灣的鹽水入侵長度進(jìn)行模擬預(yù)測，研究發(fā)現(xiàn)海平面上升1 m將使得10 ppt等鹽線向上游推進(jìn)10 km，這將改變現(xiàn)有的水質(zhì)梯度，威脅當(dāng)?shù)仫嬘盟?yīng)。Hong等[7]利用EFDC模型研究美國東海岸切薩皮克灣的河口鹽度與物質(zhì)輸運(yùn)過程對海平面上升的響應(yīng)，模擬結(jié)果顯示，海平面上升將使水體分層加強(qiáng)，同時引起鹽度上升，鹽水入侵長度增加，并且影響程度與淡水流量有明顯關(guān)系。孫志林等[8]利用FVCOM模型對錢塘江河口進(jìn)行模擬，以2012年為基準(zhǔn)年，模擬結(jié)果表明海平面上升1 m的情況下，大潮期0.45 ppt等鹽線向上游推進(jìn)了6 km，小潮期推進(jìn)了9.8 km。陳維等[9]預(yù)測海平面上升1 m時，1月長江河口的南支受北支倒灌過來的咸水的影響將明顯增大，南支水域平均鹽度均超過0.45 ppt。

海平面上升對珠江河口鹽水入侵的影響也非常顯著。何用等[10]利用磨刀門咸潮整體物理模型試驗(yàn)，對海平面上升引起的咸潮運(yùn)動變化特點(diǎn)進(jìn)行分析，發(fā)現(xiàn)海平面上升使磨刀門河道沿程鹽度均增大；以2012年的海平面為基準(zhǔn)，每升高10 cm最大咸界的上移距離大致為6.8 km。關(guān)帥等[11]利用FVCOM模型預(yù)測得海平面平均每上升1 cm，磨刀門水道鹽水入侵長度將增加0.85 km。

前人對珠江河口磨刀門鹽水入侵的研究報(bào)告較多[10]，而對東四口門的研究相對較少，且東四口門的鹽水入侵問題將影響到廣州、東莞、佛山、中山等重要城市超過3 000萬居民的生產(chǎn)生活。鑒于此，本文采用三維水動力模型EFDC(Environmental Fluid Dynamic Code)模擬珠江河口區(qū)域的水動力過程，以2007年為基準(zhǔn)年，通過設(shè)定多種海平面上升幅度的模型計(jì)算模式，對珠江河口的鹽水入侵情況進(jìn)行模擬預(yù)測，分析珠江河口東四口門(虎門、蕉門、洪奇瀝、橫門)鹽水入侵受海平面上升的響應(yīng)，為相關(guān)部門應(yīng)對海平面上升背景下該區(qū)域淡水資源管理、利用提供有效參考。

1 研究方法

1.1 研究區(qū)域

珠江河口是位于中國南海北部的一個大型河口，擁有范圍廣闊且復(fù)雜的河網(wǎng)，共有8個入?？陂T，分別是：虎門、蕉門、洪奇瀝、橫門、磨刀門、雞啼門、虎跳門和崖門，形成“三江匯流，八口入?！钡乃堤卣鱗11]，而東四口門是指八大口門中向東注入伶仃洋的4個口門[12]。珠江河口位于亞熱帶季風(fēng)氣候區(qū)，夏季高溫多雨是洪季，冬季溫和少雨是枯季。其承接西江、北江、東江三江的徑流量，多年平均總徑流量為281.1億m3[13]。

為了得到真實(shí)的水動力過程，本文將模擬區(qū)域由珠江河口擴(kuò)展至其鄰近海域，模擬網(wǎng)格見圖1。

a)模型網(wǎng)格

b)珠江河口圖1 珠江河口及鄰近海域的模擬網(wǎng)格

珠江河口的東四口門區(qū)域見圖2，該區(qū)域河網(wǎng)復(fù)雜，有多條支流，其中R1—R12是東四口門的12條主要支流，圖上還標(biāo)注了東四口門的4個自來水廠的位置，即新塘水廠、石溪水廠、沙灣水廠和東涌水廠。

圖2 東四口門及其主要河網(wǎng)

1.2 研究方法

本文采用的三維水動力模型EFDC是由美國威廉瑪麗大學(xué)維吉尼亞海洋科學(xué)研究所John Hamrick等人開發(fā)的三維地表水水質(zhì)數(shù)學(xué)模型[14]，可實(shí)現(xiàn)河流、湖泊、河口和海洋等水體的水動力和水質(zhì)模擬，是一個多參數(shù)有限差分模型，已經(jīng)在很多河口海岸地區(qū)有過成功應(yīng)用[15-17]。

模型利用曲線正交網(wǎng)格模擬珠江河口及其鄰近海域，垂向上有10個σ層。模型以珠江河口2007年的自然條件為基準(zhǔn)，包括地形、風(fēng)場、徑流量等。潮汐諧波數(shù)據(jù)由美國俄勒岡州立大學(xué)潮汐預(yù)測軟件(OTPS)算得，并用于模型開邊界的驅(qū)動。開邊界的溫度和鹽度數(shù)據(jù)取自世界海洋數(shù)據(jù)庫，風(fēng)力數(shù)據(jù)來源于香港機(jī)場觀測站(圖1b)采集的2007年風(fēng)場數(shù)據(jù)。河流徑流量數(shù)據(jù)采用珠江三條主要支流的上游各站點(diǎn)2007年的徑流量觀測數(shù)據(jù)，即石角站(北江)、高要站(西江)和博羅站(東江)。各個淡水源的流量分配按照珠江流域水資源保護(hù)局公布的比例進(jìn)行，見表1。

近幾十年來，全球海平面一直呈現(xiàn)上升趨勢，并且海平面上升的速度不斷增加。根據(jù)美國Climate Change Science Program(CCSP) 2009年工作報(bào)告指

表1 珠江河口各口門徑流量分配比例

出，到2100年海平面可能升高1.0～1.1 m[18]。Intergovernmental Panel on Climate Change(IPCC) 預(yù)測到21世紀(jì)末海平面將上升0.82 m[19]。從規(guī)避風(fēng)險(xiǎn)的角度來說，對未來海平面上升估計(jì)的結(jié)果應(yīng)該更偏重于預(yù)測上升值的99分位[20]，即區(qū)間內(nèi)至少有99%的概率包含真值。基于典型濃度路徑RCP8.5，到2100年全球平均海平面上升的99分位是0.4～1.8 m[21]。故本文以2007年為基準(zhǔn)，設(shè)定5種模式來研究珠江河口全年對不同海平面上升幅度的響應(yīng)，模式命名分別為Base、P30、P50、P100、P150，分別代表基準(zhǔn)情況，海平面上升30、50、100、150 cm的情況。各模式之間除了海平面高度不同以外，其余設(shè)置均相同，各模式從河道上游輸入的徑流量都按照2007年珠江河口的實(shí)際徑流量進(jìn)行。

1.3 模型驗(yàn)證

將珠江河口2007年的模型結(jié)果與實(shí)際觀測數(shù)據(jù)進(jìn)行比較，驗(yàn)證模型的可靠性。

1.3.1水位驗(yàn)證

取珠江河口區(qū)域的8個站點(diǎn)(三灶、香港、MO1、九洲港、內(nèi)伶仃、桂山島、燈籠山、黃埔，圖1b)作為水位驗(yàn)證點(diǎn)，所使用的水位驗(yàn)證實(shí)測數(shù)據(jù)均為2007年在各站點(diǎn)的實(shí)際觀測數(shù)據(jù)。本文利用Wilmott[22]提出的Skill值來反映模型模擬數(shù)據(jù)與實(shí)測數(shù)據(jù)之間的匹配程度。Skill按下式計(jì)算：

若模型模擬結(jié)果與實(shí)際觀測吻合越好，Skill值越接近于1，反之趨近于0。本文所選的8個驗(yàn)證點(diǎn)的Skill值見表2。結(jié)果表明模型的Skill值與1相近，模型模擬結(jié)果與實(shí)際觀測吻合程度較好。

表2 所選驗(yàn)證點(diǎn)的水位驗(yàn)證數(shù)據(jù)

各站點(diǎn)2007年12月的水位實(shí)測值和模擬的水位值隨時間的變化過程曲線見圖3。各個站點(diǎn)的水位實(shí)測值和模擬值都非常接近，模型模擬的水位與實(shí)際水位能較好吻合。

a)三灶 b)香港

c)MO1 d)九洲港

e)桂山島 f)燈籠山

g)黃埔 h)內(nèi)伶仃

圖3 各個站點(diǎn)水位驗(yàn)證曲線

1.3.2潮流驗(yàn)證

以圖1b中的MO1站點(diǎn)作為潮流驗(yàn)證點(diǎn)，把驗(yàn)證點(diǎn)處的2007年7月第1周的模擬數(shù)據(jù)與2007年該站點(diǎn)實(shí)測潮流數(shù)據(jù)進(jìn)行對比。驗(yàn)證點(diǎn)水深8.3、9.3 m處的潮流矢量變化情況見圖4。圖中可以看出，8.3、9.3 m處的實(shí)測潮流與模擬潮流的流速與流向都吻合良好。

a) 水深為8.3 m

b)水深為9.3 m

圖4 MO1站點(diǎn)的潮流驗(yàn)證

1.3.3鹽度驗(yàn)證

以圖1b中的MO23站點(diǎn)作為潮流驗(yàn)證點(diǎn)，驗(yàn)證該處在珠江典型枯水期模型模擬的鹽度變化情況與實(shí)際觀測結(jié)果是否相符。具體驗(yàn)證的時間為2007年第359.5～361.5天，驗(yàn)證點(diǎn)表層、中層和底層鹽度變化見圖5，驗(yàn)證結(jié)果反映模型模擬的鹽度在表層、中層和底層都與實(shí)際觀測的鹽度非常吻合，模型模擬效果良好，滿足使用要求。

a) 表層圖5 MO23站點(diǎn)的鹽度驗(yàn)證

b) 中層

c) 底層

續(xù)圖5MO23站點(diǎn)的鹽度驗(yàn)證

2 結(jié)果

2.1 鹽度變化

珠江河口各個口門的上游就是珠江三角洲城市群，若出現(xiàn)鹽水入侵將對周邊區(qū)域居民的生產(chǎn)生活造成嚴(yán)重影響。模型模擬的結(jié)果包含整個模擬區(qū)域全年的4D鹽度數(shù)據(jù)。通過對鹽度進(jìn)行篩選，對各口門1年內(nèi)鹽度超過1 ppt的天數(shù)進(jìn)行統(tǒng)計(jì)，得到了各個口門在不同海平面上升幅度情況下，鹽度超過1 ppt的總天數(shù)，見圖6。海平面上升對珠江河口各個口門的鹽度有明顯的影響，海平面上升幅度越大，鹽度超過1 ppt的天數(shù)越多。海平面上升1 m時，各個口門的高鹽天數(shù)相比基準(zhǔn)情況增多的天數(shù)統(tǒng)計(jì)見表3。

圖6 天數(shù)統(tǒng)計(jì)

口門增加天數(shù)/d虎門63蕉門103洪奇瀝60橫門39磨刀門110雞啼門84虎跳門155崖門153

圖6和表3中的結(jié)果表明，在相同海平面上升幅度的影響下，不同口門的響應(yīng)情況是不同的，其中虎跳門和崖門對海平面上升最敏感，海平面上升1 m時高鹽天數(shù)相比基準(zhǔn)情況分別增加155、153 d。通過研究各口門鹽度超過1ppt情況的出現(xiàn)時間，發(fā)現(xiàn)在海平面上升的情況下，各口門高鹽現(xiàn)象仍是多發(fā)生于春、冬兩季，也是珠江的枯水期，這與基準(zhǔn)情況相似，但高鹽現(xiàn)象的持續(xù)時間會隨海平面上升幅度的增加而延長。

當(dāng)自來水廠取水口處水體的鹽度超過自來水取用標(biāo)準(zhǔn)(0.5 ppt)時，自來水廠將可能采取停止自來水供應(yīng)的措施。珠江上游區(qū)域4個自來水廠，水廠位置見圖2，全年最高鹽度和超標(biāo)天數(shù)的統(tǒng)計(jì)見表4。將海平面上升1 m的模擬結(jié)果與基準(zhǔn)情況進(jìn)行對比可以看出，各個水廠鹽度受海平面上升影響非常明顯，鹽度超標(biāo)現(xiàn)象開始出現(xiàn)，尤其是東涌水廠的最高鹽度將增加7.84 ppt，同時鹽度超標(biāo)天數(shù)增加了114 d。

表4 天數(shù)統(tǒng)計(jì)

2.2 鹽水入侵長度變化

珠江河口各口門上游區(qū)域河流分叉較多，單以某一條支流的鹽水入侵長度來代表一個口門的鹽水入侵并不充分，故將各個口門的主要支流(圖2)都作為研究對象進(jìn)行分析。模擬的珠江河口在不同海平面上升幅度情況下，2007年1月各主要支流平均的0.5 ppt等鹽線位置見圖7。海平面上升幅度越大，0.5 ppt等鹽線向上游推進(jìn)越多，鹽水會入侵到越上游區(qū)域，加劇鹽水入侵。

圖7 平均0.5 ppt 等鹽線所在位置

海平面上升對珠江河口的鹽水入侵有明顯的強(qiáng)化作用。基準(zhǔn)情況下，洪季各口門都不出現(xiàn)鹽水入侵，但當(dāng)海平面上升150 cm時，即使在洪季，各口門也會出現(xiàn)不同程度的鹽水入侵問題。把圖2中標(biāo)定的各口門位置作為計(jì)算鹽水入侵長度的起算點(diǎn)，以各支流0.5 ppt等鹽線沿河道至對應(yīng)口門的距離作為鹽水入侵長度。海平面上升100、150 cm時，東四口門各支流鹽水入侵長度相比基準(zhǔn)情況向上游的最大推進(jìn)距離見表5。由于蕉門(R7、R8)、洪奇瀝(R9、R10)和橫門(R11、R12)各處鹽水都未入侵到各自的河流分叉點(diǎn)(圖2中A、B、C為分叉點(diǎn))，故這3個口門各自的支流的入侵情況相同。表 5的結(jié)果可以看出海平面上升1 m時，各支流的0.5 ppt等鹽線相較于基準(zhǔn)情況都向上游推進(jìn)了較大距離，其中R2的推進(jìn)距離最大，達(dá)21.371 km；蕉門最大推進(jìn)9.644 km；洪奇瀝最大推進(jìn)9.753 km，橫門最大推進(jìn)4.819 km。

表5 0.5 ppt等鹽線相比基準(zhǔn)情況最大向上游推進(jìn)的距離 km

3 討論

受Yang等[23]和Chen等[24]將海平面上升幅度與鹽度變化通過數(shù)據(jù)回歸的方式得到擬合度較好的直線或曲線的啟發(fā)，并考慮到本文模型設(shè)定的幾個計(jì)算模式之間，僅有海平面上升幅度一個變量，便將各條支流鹽水入侵長度進(jìn)行月平均，與海平面上升幅度進(jìn)行線性回歸，得到海平面上升幅度與各支流鹽水入侵長度的擬合直線。選擇入侵情況最嚴(yán)重的1月份作為線性回歸分析對象，得到的回歸直線(圖8)可以用于預(yù)測未來海平面上升某一特定幅度時，珠江河口各口門上游區(qū)域鹽水入侵的距離。每條回歸直線的誤差范數(shù)見表6。誤差范數(shù)都很小，則回歸直線可靠，說明珠江河口各支流由于海平面上升幅度變化引起的鹽水入侵長度變化與海平面上升幅度之間存在著較好的線性關(guān)系。得到的線性回歸關(guān)系可用于預(yù)測特定海平面上升幅度下某支流的鹽水入侵長度，指導(dǎo)工程實(shí)踐，如自來水廠遷址、選址等。

表6 R1-R12線性回歸直線誤差范數(shù)

a)虎門R1 b)虎門R2 c)虎門R3

d)虎門R4 e)虎門R5 f)虎門R6

g)蕉門R7-R8 h)洪奇瀝R9-R10 i)橫門R11-R12

圖8 海平面上升幅度與鹽水入侵長度的回歸直線

4 結(jié)論

海平面上升會加劇珠江河口東四口門的鹽水入侵問題。在海平面上升1 m的情況下，東四口門高鹽天數(shù)(鹽度≥1 ppt)都有明顯增加，其中蕉門增加最多，達(dá)3個多月。而建在東四口門上游的自來水廠也會出現(xiàn)由于海平面上升引起的取水口處鹽度超標(biāo)的問題，其中新塘水廠和東涌水廠受到的影響最大，海平面上升1 m情況下，全年鹽度超標(biāo)天數(shù)將增加3個多月，所以在未來各個水廠必須要考慮遷址。

模型模擬結(jié)果還表明海平面上升會使東四口門的鹽水入侵長度增加，其中虎門各支流鹽水入侵長度增加最多，海平面上升150 cm情況下，0.5 ppt等鹽線最多將向上游推進(jìn)27 km。通過對珠江河口東四口門主要支流在一月份的平均鹽水入侵長度與海平面上升幅度進(jìn)行線性回歸，得到了可靠的回歸直線。借助回歸直線，可以預(yù)測未來海平面上升特定幅度時，珠江河口各支流的鹽水入侵長度，可以為未來水廠遷址提供依據(jù)。