西北江三角洲河床變化對咸潮上溯影響的數值模擬研究

http：//www.renminzhujiang.cn

DOI：10.3969/j.issn.1001-9235.2024.03.010

靳高陽，朱三華，曹一豪.西北江三角洲河床變化對咸潮上溯影響的數值模擬研究［J］.人民珠江，2024，45（3）：89-98.

摘"要：為研究西北江三角洲河床變化及其與咸潮上溯關系，收集整理以往研究成果資料與水文氣象、河道地形數據的基礎上，建立一維非恒定流含氯度數學模型，水動力和含氯度驗證結果均滿足相關規范要求，可以用于計算分析河床變化對咸潮上溯的影響。研究結果表明：①1999—2016年，西北江三角洲網河區河床處于不均勻下切狀態，河道容積持續增大，導致河口漲落潮阻力減小，高鹽水團沿河口上溯更為簡單，咸潮上溯情勢日益嚴峻；②各口門站的漲潮量增加，咸界上移取淡機率降低，思賢滘不同來流條件下，磨刀門水道咸界平均上移3 735 m，平崗泵站的取淡機率平均減少5.33%，沙灣水道咸界平均上移2 369 m，沙灣第一水廠取淡機率平均減少5.14%，嚴重影響河口地區的生態系統以及正常取水需求。研究結果可以為西北江三角洲及河口地區的供水安全保障措施研究及工程設計提供基礎支撐。

關鍵詞：河床演變；含氯度；數學模擬；咸潮上溯；咸界變化

中圖分類號：TV8""文獻標識碼：A""文章編號：1001-9235（2024）03-0089-10

Numerical Simulation Study on Influence of Riverbed Changes on Saltwater Intrusion in Northwest River Delta

JIN Gaoyang1，ZHU Sanhua1，CAO Yihao2

（1.China Water Resources Pearl River Planning，Surveying amp; Designing Co.，Ltd.，Guangzhou 510610，China;

2.School of Civil Engineering，Sun Yat-sen University，Guangzhou 510275，China）

Abstract：To study the changes in the riverbed of the Northwest River Delta and their relationship with saltwater intrusion，this paper builds a one-dimensional unsteady flow chlorine concentration mathematical model based on the collection and organization of previous research results，hydrological，meteorological，and river terrain data.The hydrodynamic and chlorine concentration verification results meet the relevant regulatory requirements and can be applied to calculate and analyze the impact of riverbed changes on saltwater intrusion.The research results are as follows.① From 1999 to 2016， the riverbed of the Northwest River Delta network was in an uneven downward cutting state，and the volume of the river continued to increase，resulting in a decrease in the resistance of the estuary to rising and falling tides.The upward tracing of high saltwater masses along the estuary was simpler，and the upward tracing of salt tides became increasingly severe.② The rising tide at each Koumen station increases，and the probability of removing freshwater from the saltwater boundary decreases.Under different inflow conditions in Sixianjiao，the average upward movement of the saltwater boundary at Modaomen Waterway is 3 735 m，and the average downward movement of the freshwater from Pinggang Pumping Station is 5.33%.The average upward movement of the saltwater boundary at Shawan Waterway is 2 369 m，and the average downward movement of the freshwater from Shawan No.1 Water Plant is 5.14%.These seriously affect the ecosystem and normal water intake demand in the estuarine area.The research results can provide basic support for the research and engineering design of water supply safety measures in the Northwest River Delta and estuarine areas.

Keywords：riverbed evolution;chlorine content;mathematical simulation;saltwater intrusion;saltwater changes

咸潮上溯是指當河口淡水流量相對較小時，海洋高鹽水團在太陽和月球對地表水的引力作用下，隨潮汐周期性沿河口上溯造成河道水體變咸的天然水文現象^［1］。一般來說，咸潮上溯會對河口地區的生態系統以及正常取水需求產生較大的影響。咸潮上溯作為世界河口的普遍性現象，是國內外學者的重要研究對象：Savenije^［2］通過對大量實測資料的分析，建立了咸潮上溯曲線與河口的幾何相關關系，并指出幾何參量是影響河口咸潮上溯數學模型精度的關鍵因子；而Gualbert等^［3］則采用MOCDENS 3D 模型模擬分析了三維空間條件下河口地區的咸潮上溯規律。國內學者對咸潮上溯現象的研究則主要集中在長江河口和珠江河口，主要研究方法可分為對實測資料的數據統計分析以及數值模擬分析兩大類^［4-5］。有研究通過多年的資料指出，河口咸潮上溯主要受上游徑流量、三角洲河床地形演變、溫度以及海洋潮汐動力等控制^［6-9］；在此基礎上，基于實測的現場資料統計分析了珠江河口咸潮上溯的現狀和成因，并明確給出影響咸潮的關鍵因子^［10-11］。近年來，在統計分析的基礎上，咸潮上溯數值模型逐漸建立：黃銳貞等^［7］通過建立基于BP神經網絡的含氯度變化模擬模型，分析了上游徑流量變化對含氯度的影響及不同來水頻率和潮位設計方案下的含氯度變化；劉祖發等^［12］利用FVCOM模型對徑流量變化引起的磨刀門水道咸潮上溯現象進行數值模擬，厘清了上游河道水位上升后水動力的變化規律；齊慶輝等^［13］則是基于MIKE建立了覆蓋珠江八大口門區的高精度三維水流鹽度數學模型，模擬出珠江口地區三維水動力格局與鹽度分布。

珠江三角洲地區作為中國重要的經濟中心區域，其河口水源是維系周邊城市工農業生產和人民生活的重要保障。但近年來，隨著用水量的大幅提高、枯季干旱、三角洲河道地形演變、枯季西北江分流比變化等多重因素影響，咸潮上溯呈加劇趨勢，咸界上移^［14］，甚至在2022年出現了磨刀門水道取水口自12月4日起全線24 h含氯度超標，其中珠海的平崗泵站連續10 d含氯度超標，連續半月無法取水的情況，人民的生產生活受到了嚴重影響。因此，厘清咸潮上溯成因及其作用機理就顯得尤為重要。但由于珠江三角洲咸潮上溯因子錯綜復雜：上游徑流量、河道下切度、外海的波浪與風等都會引起咸潮上溯距離的變化，且珠江三角洲口門眾多，因此，本文在收集整理以往研究成果與現有水文氣象、河道地形的基礎上，利用實測地形與以往水文同步測驗資料建立一維非恒定流含氯度數學模型，研究西北江三角洲河床變化對咸潮上溯的影響，為西北江三角洲及河口地區的供水安全保障工程深入研究及設計提供基礎支撐。

1"含氯度數學模型的建立與驗證

1.1"一維非恒定流含氯度模型

珠江河網一維非恒定流數學模型采用半隱式有限元法進行計算。其水流運動的基本方程由一維水流方程組及汊點連接方程組兩部分組成，具體方程如下。

1.1.1"河道流動的圣維南方程組

連續方程：

式中"x、t——空間、時間坐標；A——河道過水斷面面積；Q——斷面流量；q——均勻旁側入流；Q_c——集中旁側入流；δ——Dirac函數；Z——水位水量；K——流量模數，由謝才公式計算^［15］。

由式（1）、（2），可得含氯度對流－擴散運動方程如下：

式中"s——水體含氯度；Q_s——氯化物量；D——紊動擴散系數。

含氯度的對流擴散系數的確定是保證含氯度數學模型精度的關鍵。通常情況下，擴散系數可按紊流比擬概念計算，根據實測資料分析表明，采用常用的紊流模型計算的擴散系數比資料顯示的擴散系數至少小一個量級。因此，本文依據量綱分析法，采用如下更為準確的紊流擴散模型^［15］：

D=（ν_ts+ν_a）/Sch（4）

式中"ν_ts——Smagorinsky紊流渦黏系數；ν_a——潮流與風浪引起的附加紊流擴散系數；Sch——Schimit數。

ν_a表達式如下：

式中"H——潮差；g——重力加速度；C_w——風浪影響因子，由實測含氯度計算標定。

1.1.2"河道單元的離散^［15］

1.1.2.1"連續方程的離散

連續方程的對流方程形式見式（6）：

式中，N、W分別為插值函數與權函數。經有限元空間半離散后的方程為常微分方程，可采用多種方法求解。注意到過水斷面面積A是水位Z的函數，因此，連續方程的單元離散方程可轉化為求解水位的方程。

1.1.2.2"運動方程的單元離散

運動方程的對流方程見式（8）：

1.1.3"河網汊點連接方程^［15］

1.1.3.1"質量守恒關系

進出每一汊點必須與該汊點蓄水量的增減相平衡，即節點的質量守恒方程：

式中"Ω、Z——汊點的蓄水量、水位；A_c——汊點的蓄水面積（匯合區面積）；Q_i、Z——通過I河道斷面進入該汊點的流量、汊點水位。

1.1.3.2"水位銜接關系

節點一般可概化成一個幾何點，出入各汊節點的水位平緩，不存在水位突變情況，則各節點相連汊道的水位應相等，等于該點的平均水位，即Z=Z_i。

綜上，通過對連續方程和運動方程中流量和水位的求解，可以得到一維非恒定流含氯度數學模型中的水體含氯度。

1.2"模型范圍與斷面布置

一維含氯度數學模型的計算范圍，上起西江的梧州站、北江的飛來峽站、流溪河的老鴉崗站、增江的麒麟咀站、東江的博羅站，下至伶仃洋舢舨洲、蕉門延伸段出口、洪奇門和橫門匯合延伸段出口、橫門南汊出口、磨刀門馬騮洲和橫洲口、雞啼門出海口和黃茅海，河道總長度2 091 km，計算網絡包括266個汊口，403個河段，2 388個河道斷面，33個邊界。模型主要采用1999—2016 年珠江流域西北江河道地形實測數據，由于部分外海區域由于缺乏最新的地形資料，本次模擬采用1999、2016年2套實測地形資料。西北江三角洲一維含氯度數學模型范圍及斷面位置示意見圖1。

1.3"數學模型驗證

本次含氯度計算，主要分析河道地形變化導致珠江三角洲枯水期含氯度變化情況，1999年地形采用“2005·1”枯季實測的水文資料對模型進行流量率定，采用“2005·1”實測含氯度資料對含氯度參數進行率定。2016年地形采用“2019·12”枯季實測的水文資料對模型進行流量率定，采用“2020·1”實測含氯度資料對含氯度參數進行率定。

1.3.1"水動力驗證

本次模型驗證過程中，2016年地形河道枯水流量采用“2019·12”枯水進行率定（12月27日14：00至1月5日14：00），共217 h。1999年地形枯水河道流量采用“2005·1”枯水進行率定（1月18日19：00至2月5日7：00），共421 h。由于本模型驗證圖較多，本文只節選2019年12月枯季條件下各站點流量驗證過程（圖2、3）；珠江三角洲主要站點特征潮位率定結果見表1。

由圖2和表1，可以得到：主要站點計算潮位特征值與實測潮位相比，誤差均小于0.10 m，計算潮位過程與實測潮位過程吻合，相位誤差小于1 h；主要站點斷面計算最大漲落潮流量與實測最大漲落潮流量誤差一般小于10%，基本滿足JTJ/T 233—98《海岸與河口潮流泥沙模擬技術規程》的規范要求，說明模型采用的斷面和網格剖分、糙率等參數基本合理。

1.3.2"含氯度驗證

本次含氯度驗證采用2020年1月17日18時至2月1日17時（共360 h）實測含氯度資料對模型的擴散參數進行率定，思賢滘的平均流量2 400 m³/s。2020年1—2月枯季含氯度驗證過程見圖4、5，由圖可知，含氯度計算與實測數據吻合良好，說明采用擴散參數合理，含氯度模型可靠，可以用于含氯度規律的研究。

2"咸潮上溯規律分析

近20 a來，由于珠江三角洲及港澳地區快速發展和大規模建設，導致珠江三角洲地區出現大規模的采砂活動，疊加珠江三角洲高等級航道建設，造成三角洲河床大規模下切，漲落潮阻力減小，漲潮空間增大，咸潮上溯距離增加。但由于河床下切不均衡，導致西、北江水道的分流比發生變化，各口門壓咸流量相應發生變化，進而造成咸潮上溯規律錯綜復雜。本文采用1999年地形和2016年地形，計算分析河床變化及其與咸潮上溯關系。

2.1"不同年代地形變化規律

1999—2016年，珠江三角洲河道容積增加共70 129萬m³，其中西江片河道容積增加40 504萬m³，河床平均下切2.01 m，北江片河道容積增加29 625萬m³，河床平均下切1.06 m。經分析，大范圍大幅度的河床下切，主要與航道疏浚、河道采砂等人類活動有關。

2.2"不同年代地形水動力變化規律

本文采用一維非恒定流數學模型，上邊界按思賢滘流量2 500 m³/s（三水+馬口）控制，東江博羅站、增江麒麟咀站、流溪河老鴉崗站和潭江石咀站采用2019年12月枯季實測的來流，模型下邊界采用2019年12月枯季實測的水位，在珠江三角洲1999年地形和2016年地形條件下，計算得到的低低潮位水面線沿程變化見圖6、7。

由圖6、7可知，1999—2016年，西北江三角洲主要河床一直處于下切狀態，河道容積持續增大。北江河床（北江干流-順德水道-沙灣水道）平均下切1.51 m，西江河床（西江干流-西海水道-磨刀門水道）平均下切2.44 m，西北江三角洲河床落潮順暢，低潮位平均降低0.09 m。多年來西北江三角洲河床的持續下切導致河口漲落潮阻力減小，高鹽水團沿河口上溯更為簡單，咸潮上溯情勢日益嚴峻。

落潮平均流速影響著河口區干流水體鹽度的變化過程，是分析咸界位置和取淡幾率的重要因子。在珠江三角洲1999年地形和2016年地形條件下，利用一維非恒定流模型計算得到的西北江干流平均流速沿程變化見圖8、9。

圖8、9表明：對西江干流來說，不同年代地形條件下西江干流落潮平均流速呈略微減小趨勢，但整體變化不大：在1999年地形條件下，西江干流落潮平均流速約為0.32 m/s，而在2016年地形條件下西江干流落潮平均流速約為0.31 m/s，落潮平均流速整體減少0.01 m/s；不同于西江干流，北江干流不同年代地形落潮平均流速稍有加大，1999年地形北江干流落潮平均流速為0.34 m/s，2016年地形北江干流落潮平均流速為0.36 m/s，落潮平均流速增加0.02 m/s。

2.3"不同年代地形咸界變化規律

本文采用一維非恒定流含氯度模型，模擬西北江三角洲不同年代地形、思賢滘不同流量情況下咸界變化規律，計算統計主要口門站漲潮量變化、主要水廠超標時間以及主要河道咸界上溯距離變化見表2—4。

從表2可知，西北江三角洲口門站的漲潮量，隨著上游徑流量增加，徑流動力增強，潮流動力相對減弱，漲潮量減少；且2016年地形相比1999年地形，由于三角洲河床下切，槽蓄容積增加，口門站的漲潮量增加，燈籠山漲潮量增加16.65%～19.40%，三沙口漲潮量增加3.10%～18.93%；而由表3可知，由于口門的漲潮量增加，平崗泵站的超標時間平均增加2.78 h，沙灣第一水廠超標時間平均增加2.38 h，平崗泵站的取淡機率平均減少5.33%，沙灣第一水廠取淡機率平均減少5.14%，不斷影響著城市居民的正常取水需求。表4則明確指出了一維非恒定流模型下2016年地形條件相比1999的地形咸界變化情況：由于河道地形下切，口門漲潮量增加，咸界上移，磨刀門水道咸界平均上移3 735 m，沙灣水道咸界平均上移2 369 m，咸潮上溯情勢日漸嚴峻。

綜上所述，從1999—2016年，西北江三角洲由于各種人類活動的影響，導致河床的持續下切，河道槽蓄容積增加，各個口門站的漲潮量增加，導致河道上游水體整體變咸，咸界不斷上移，主要水廠的取淡幾率減小、水質超標時間增加，嚴重影響河口地區的生態系統以及正常取水需求。不同年代地形條件下思賢滘不同來流時咸界變化見圖10、11。

3"結論

a）本文在收集整理以往研究成果與現有水文氣象、河道地形的基礎上，利用實測地形與以往水文同步測驗資料建立一維非恒定流含氯度數學模型，水動力和含氯度驗證結果均滿足相關規范要求，可以用于計算分析河床變化及其與咸潮上溯關系。

b）不同年代地形變化導致枯季水動力變化規律：由于河道下切，西北江三角洲河道落潮順暢，1999—2016年，低潮位平均降低0.09 m；且不同年代地形落潮流速變化較小，落潮平均流速變化值在0.02 m/s以內。

c）不同年代地形上游不同來流時，隨著上游徑流量增加，徑流動力增強，潮流動力相對減弱，口門站的漲潮量減少。

d）從1999—2016年，由于河床下切，水深增加，漲潮阻力減少，各口門站的漲潮量增加，咸界上移取淡機率降低，思賢滘不同來流條件下，磨刀門水道咸界平均上移3 735 m，平崗泵站的取淡機率平均減少5.33%，沙灣水道咸界平均上移2 369 m，沙灣第一水廠取淡機率平均減少5.14%，嚴重影響河口地區的生態系統以及正常取水需求。研究成果可以為西北江三角洲及河口地區的供水安全保障工程深入研究及工程設計提供基礎支撐。

參考文獻：

［1］趙高磊，張志廣，梁瑞峰，等.北門江河口徑流變化對河口段水環境的影響［J］.人民長江，2018，49（15）：36-40，95.

［2］SAVENIJE H.Rapid assessment technique for salt intrusion in alluvial estuaries［R］.Delft：Hydraulics and Environmental Engineering，1992．

［3］GUALBERT H P.Oude Essink，Saltwater intrusion in 3D large-scale aquifers：a dutch case，Physics and Chemistry of the Earth［J］.Oceans and Atmosphere，2001，26（4）：337-344.

［4］LI L，ZHU J R，WU H，et al.A numerical study on water diversion ratio of the Changjiang （Yangtze） estuary in dry season［J］.Chinese Journal of Oceanology and Limnology，2010，28（3）：700-712.

［5］韓乃斌.長江口南支河段氯度變化分析［J］.水利水運科學研究，1983（1）：74-81.

［6］LIAN J J，WEI H，CHAO M，et al.Guarantee rate of freshwater in a river mouth intruded by saltwater with respect to the joint impact of runoff and tide［J］.Journal of Hydroinformatics，2015，17 （6）：917-929.

［7］黃銳貞，陳曉宏，何艷虎，等.上游來水與下游潮位對河口含氯度的影響：以磨刀門河口為例［J］.中國農村水利水電，2017（4）：101-105.

［8］沈煥庭，茅志昌，朱建榮.長江河口鹽水入侵［M］.北京：海洋出版社，2003.

［9］朱建榮，傅利輝，吳輝.風應力和科氏力對長江河口沒冒沙淡水帶的影響［J］.華東師范大學學報（自然科學版），2008，39（6）：1-8，39.

［10］宋曉飛，石榮貴，孫羚晏，等.珠江口磨刀門鹽水入侵的現狀與成因分析［J］.海洋通報，2014，33（1）：7-15.

［11］孔蘭，陳曉宏.珠江口咸潮影響因素分析［J］.水資源保護，2015，31（6）：94-97，134.

［12］劉祖發，張泳華，陳曉越，等.數值分析徑流變化對磨刀門咸潮上溯的影響［J］.亞熱帶資源與環境學報，2017，12（4）：30-38.

［13］齊慶輝，諸裕良，謝寧寧.珠江河口航道整治工程對咸潮上溯的影響研究［J］.水運工程，2015（10）：138-143.

［14］韓志遠，田向平，劉峰.珠江磨刀門水道咸潮上溯加劇的原因［J］.海洋學研究，2010，28（2）：52-59.

［15］朱三華，沈漢堃，金生，等.珠江三角洲含氯度模型研究［J］.人民珠江，2010，31（S1）：9-13.

（責任編輯：高天揚）

基金項目：國家自然科學基金（U1911204、51861125203）;國家重點研發計劃項目（2021YFC3001000）

收稿日期：2023-07-17

作者簡介：靳高陽（1984—），男，高級工程師，主要從事水利規劃與設計工作。E-mail：12584298@qq.com