基于河口三維斜壓數(shù)學模型的珠江河口咸淡水比例研究

黃廣靈，邱 靜，譚 超，黃本勝

(1.廣東省水利水電科學研究院，廣州 510635；2.廣東省流域水環(huán)境治理與水生態(tài)修復重點實驗室，廣州 510635；3.河口水利技術(shù)國家地方聯(lián)合工程實驗室，廣州 510635)

1 概述

世界很多發(fā)達城市均位于河口地區(qū)[1]，河口區(qū)復雜的水動力交匯過程給河口水資源利用帶來與內(nèi)陸河流不一樣的問題，如咸潮對淡水資源的影響等，因此，合理利用及管理河口淡水資源是解決工業(yè)及城市居民生活用水的關(guān)鍵問題。珠江河口三角洲網(wǎng)河是粵港澳大灣區(qū)重要水源區(qū)[2]，然珠江三角洲受潮汐作用影響，下游河口區(qū)淡水資源受到近岸高鹽陸架水入侵的“污染”[1]，致使灣區(qū)局部地區(qū)淡水資源的利用受到較大限制，特別是對水質(zhì)鹽度要求的生活、生產(chǎn)取水，每年枯季均會有數(shù)天因取水口鹽度超標而無法實現(xiàn)取水，嚴重影響大灣區(qū)特別是下游區(qū)的取水保障率。因此，從河口區(qū)可利用水資源量的角度，需要對珠三角感潮區(qū)的咸淡水(咸淡水指河口區(qū)含氯度大于250 mg/l的微咸水)進行定量計算及分析。

對于珠江河口咸淡水，以往研究主要關(guān)注咸潮的入侵規(guī)律[3-6]、影響因素[7]、預報技術(shù)[8]以及抑咸措施[9-11]咸潮活動自身的特征，對咸潮活動影響下的河口水資源利用研究目前仍較少，且也僅限于定性分析[12]，缺乏定量的計算。隨著最嚴格水資源管理制度的不斷深化，咸淡水作為國家鼓勵的非常規(guī)水資源，如何計入用水總量控制指標，均直接關(guān)系到水資源“三條紅線”控制及其考核管理，尤其是取水許可管理、取用水量指標統(tǒng)計、水資源費征收等諸多水資源管理問題。計算咸淡水的比例成為河口區(qū)咸淡水管理的核心問題，目前的咸淡水比例計算常規(guī)方法是根據(jù)水文站長序列鹽度監(jiān)測資料進行統(tǒng)計，但存在受水文站點布設(shè)制約，不能全面反映河口區(qū)不能水域的咸淡水比例，存在較大的局限性。因此，本文將采用數(shù)學模型計算，研究珠江河口年尺度下的咸潮變化情況，以此定量分析河口咸淡水資源量，以為河口區(qū)咸淡水資源的管理提供科學依據(jù)。

2 三維斜壓數(shù)學模型

本文采用SELFE[13]模式構(gòu)建了珠江網(wǎng)河—河口灣—近海陸架整體三維斜壓數(shù)學模型。

2.1 控制方程

模型基于靜壓近似及Boussinesq近似求解三維淺水方程和溫鹽輸送方程。 模型求解的主要未知變量包括：自由水面高程、流體速度、溫度和鹽度矢量。在笛卡爾坐標系統(tǒng)下，模型的主要控制方程如下：

(1)

(2)

(3)

(4)

(5)

(6)

式中：

?——哈密頓算子；

(x,y)——水平笛卡爾坐標，m；

z——垂向坐標，向上為正，m；

t——時間，s；

η(x,y,t)——自由水面高程，m；

h(x,y)——水深，m；

w——垂向流速，m/s；

f——柯氏力因子，S-1；

g——重力加速度，m/s2；

ψ(φ,λ)——潮汐勢，m；

α——有效地球彈性因子；

ρ(x,t)——水的密度，默認參考值ρ0為1 025 kg/m3；

pA(x,y,t)——自由水面大氣壓強，N/m2；

S,T——水的溫度和鹽度，psu；

υ——垂向渦動粘性系數(shù)，m2/s；

μ——水平渦動粘性系數(shù)，m2/s；

κ——溫度和鹽度的垂向渦動擴散系數(shù)，m2/s；

Fs,Fh——輸移方程中的水平擴散系數(shù)；

Q——太陽輻射吸收率，W/m2；

Cp——水的比熱，J/(kg/K)。

模型中所有方程采用半隱格式離散。連續(xù)方程及動量方程同時求解，通過底邊界條件，將連續(xù)方程(1)及動量方程(3)耦合。動量方程(3)中的平流項采用歐拉-拉格朗日(ELM)方法求解，而物質(zhì)輸移方程(4)(5)則可采用ELM或者迎風有限體積法求解。

2.2 模型范圍及網(wǎng)格

建立范圍覆蓋了珠江三角洲-河口-外海一體的咸潮計算數(shù)學模型。模型上邊界選取西江高要，北江石角，東江博羅，流溪河及潭江上游，外海下邊界最外取100 m等深線處，模型網(wǎng)格如圖1所示。模型采用全三角形的非結(jié)構(gòu)性網(wǎng)格，以精確的貼合復雜的河岸邊界，模型網(wǎng)格有10萬多個節(jié)點，17萬多個單元，網(wǎng)格單元的尺度跨度為0.01～1.2 km。模型在垂向上采用Sigma坐標，共分20層，均分水深，分層厚度隨水深而變化。模型中珠江三角洲網(wǎng)河所用地形采用2005—2010年的大范圍組合地形，河口灣及近岸海區(qū)地形采用2000—2015年海圖數(shù)字化成果，外海地區(qū)地形采用 ETOPO1全球海洋地形數(shù)據(jù)。

圖1 珠江河口全三維模型網(wǎng)格示意

2.3 模型范圍及網(wǎng)格參數(shù)設(shè)置

模型采用斜壓模式，其中外海鹽度邊界取海水鹽度34 psu，上游開邊界鹽度均設(shè)為0。水平網(wǎng)格內(nèi)部各節(jié)點的初始鹽度采用以下方法設(shè)定：海域節(jié)點采用WOA09鹽度資料進行插值，網(wǎng)河節(jié)點則直接設(shè)為0。模型的干濕判定水深選取0.01 m。湍流閉合模型中動量方程和輸移方程中的水平及垂向渦擴散系數(shù)給常數(shù)10-6。在底邊界摩擦力的處理上，本文選用二次的阻力公式，網(wǎng)格中各節(jié)點的阻力系數(shù)分河段設(shè)定，并通過模型數(shù)次率定結(jié)果進行調(diào)整。

模型計算的基準面為珠基，計算時間步長為 180 s，網(wǎng)格內(nèi)各節(jié)點初始水位均設(shè)為0 m，模型計算50 h后水位穩(wěn)定，計算約8 d后鹽度場達到穩(wěn)定。模型計算水動力邊界條件包括外海及上游，其中外海是潮汐調(diào)和常數(shù)，調(diào)和常數(shù)由中國海域潮汐預報系統(tǒng)Chinatide計算得到，共選取了珠江口海區(qū)9個主要分潮：Q1，O1，P1，K1，N2，M2，S2，K2、SA。模型上游邊界均采用流量邊界。

2.4 模型率定和驗證

模型主要選取了2001年2月的珠江三角洲河網(wǎng)區(qū)同步測量水文資料進行潮位、流速、鹽度的率定，率定參數(shù)包括底床糙率和擴散系數(shù)，糙率的設(shè)置對潮位、流速起主要影響作用，而擴散系數(shù)則主要影響鹽度模擬結(jié)果，擴散系數(shù)越大則鹽水入侵范圍越快越大。模型提前于率定數(shù)據(jù)起算，總計算時長為30 d，選擇后8 d(2001年2月7日14:00—2001年2月14日20:00)進行水位、流速及鹽度率定。為了使模型更好地模擬珠江三角洲網(wǎng)河的潮汐漲落過程，選擇了八大口門的主要測站及網(wǎng)河區(qū)的各潮位站進行潮位率定。率定站點共計33個，覆蓋整個珠江三角洲。各主要測站的計算潮位與實測潮位的對比如圖2所示。率定結(jié)果顯示，絕大部分測站的潮位誤差均在0.1 m以內(nèi)，模型能夠較好的反映真實的原型特征。

圖2 潮位率定結(jié)果示意(2001年2月)

模型共選擇網(wǎng)河區(qū)及口門處共17個測站的流速資料進行率定，主要率定結(jié)果如圖3所示，從流速的率定和驗證結(jié)果來看，大部分測站的計算流速與實測流速吻合得較好，僅有少數(shù)幾個站點的流速誤差稍大，但總的誤差不超過10%，考慮到珠江網(wǎng)河及地形的復雜性，模型的流速計算結(jié)果可以較真實地反映水流的實際情況。

圖3 流速率定結(jié)果示意(2001年2月)

圖4為模型率定計算(計算時間為2001年2月)中各口門主要測站的鹽度計算值與實測值的對比；從圖4中可以看出：鹽度的模擬精度相對于潮位及流速有所降低，但總的來說計算鹽度的誤差在合理范圍以內(nèi)，各站的鹽度平均誤差均不超過實測范圍的20%[14]，能夠較好地體現(xiàn)出鹽度的大小潮變化特征?？紤]到地形、風等因素對鹽度的變化有很大影響，而模型地形與率定驗證時段地形有所差異，同時也沒有把風場加以考慮，故計算結(jié)果存在充許范圍內(nèi)的誤差符合實際情況。

圖4 鹽度率定結(jié)果示意(2001年2月)

選擇了2005年1月的實測潮位及流速對模型進行驗證，另選擇了2005年1月、2008年1月、2009年12月和2011年1月共4組鹽度實測資料對模型的鹽度計算進行驗證。在4組驗證計算中，驗證結(jié)果均顯示模型參數(shù)選擇是合理的，各計算要素的誤差均在允許范圍之內(nèi)(見圖5)。綜上可知，利用本次建立的模型進行長時間序列咸潮入侵計算合理可行。

圖5 鹽度驗證結(jié)果示意(2005年1月)

2.5 計算邊界選擇

珠江河口咸淡水比例變化受上游徑流影響較大，受枯季流量的影響尤為顯著，因此，上游流量過程采用1960—2009年每年10月—次年3月的平均流量進行排頻，計算枯季平均流量保證率，本文珠江河口咸淡水比例計算的主要典型年選取西江下游控制站高要站+北江下游控制站石角站50%、75%、90%枯季來水保證率、東江下游控制站博羅50%、75%、90%枯季來水保證率。

采用西江高要水文站、北江石角水文站及東江博羅水文站1960—2009年共50 a的流量資料，分別將西江下游控制站高要站+北江下游控制站石角站歷年枯季平均流量(以下簡稱“高+石”)、博羅站歷年枯季平均流量以P-Ⅲ適線法進行頻率計算，得到高+石及博羅站不同枯水保證率對應(yīng)流量及典型年，在此基礎(chǔ)上選擇年內(nèi)洪枯季流量差異最大的典型年作為計算邊界，最終選取的典型年見表1所示，其余上邊界流溪河老鴉崗及潭江石咀流量典型年的選取與高+石一致。

表1 不同枯水保證率對應(yīng)流量及典型年 m3/s

3 結(jié)果分析

珠江河口水資源管理對象主要為各取水戶，其取水規(guī)模相對較穩(wěn)定，因此，其取水中咸淡水的占比可以轉(zhuǎn)化成河道咸淡水時間占比計算。珠江河口河道內(nèi)咸淡水時間比例(簡稱咸淡水比例)數(shù)學模擬計算采用方法如下：本次計算的上邊界水文條件共采用3組，分別是上游枯季來水保證率為50%、75%和90%的年份的全年逐日實測流量。計算時長為1 a，從典型年的前1 a的10月到該典型年份的9月，并提前30 d起算，以保證鹽度場達到穩(wěn)定狀態(tài)。計算結(jié)果的處理：計算結(jié)果可得到每一個點的全年逐時含氯度過程，從而可統(tǒng)計每個點全年咸淡水的比例，按全年含氯度大于250 mg/L的小時數(shù)除以全年總小時數(shù)得到結(jié)果百分率統(tǒng)計，如式(7)所示。各點咸淡水比例經(jīng)空間線性插值后在繪出其咸水比例的等值線。

(7)

式中:

Sper——計算點的咸淡水比例；

Ts——該點計算年中含氯度大于250 mg/L的小時數(shù)；

Ty——計算年的總小時數(shù)。

3.1 計算咸淡水比例與實測數(shù)據(jù)的對比

3.1.1磨刀門水道

表2為枯季平均流量保證率為50%、75%和90%的水文條件下磨刀門水道掛定角、大涌口、廣昌泵站、燈籠山站及平崗泵站的年咸淡水比例計算結(jié)果。與根據(jù)實測資料的統(tǒng)計結(jié)果對比可知，在同樣的枯季平均流量保證率下通過模型計算得到的咸淡水比例更高。如根據(jù)馬口水文站枯水平均流量枯水保證率的統(tǒng)計結(jié)果，2005年高+石的流量保證率為74%，2005年的掛定角、大涌口水閘和平崗泵站的實測年咸淡水比例分別為51%、41%和13%，而高+石枯季平均流量保證率為75%水文條件模型計算得到掛定角、大涌口水閘和平崗泵站3站咸淡水比例分別為62.0%、48.0%和9.0%。

表2 磨刀門水道監(jiān)測站點咸淡水比例計算結(jié)果與實測結(jié)果 %

模型計算結(jié)果與實測值相比存在一定的誤差，主要有下幾點原因：① 河口咸淡水混合受諸多因素影響，而模型模擬過程未能考慮所有影響因子，只考慮了其中最重要的，故計算結(jié)果存在一定誤差是合理的；② 與數(shù)據(jù)的處理方法有關(guān)，本文中所出現(xiàn)的咸淡水比例結(jié)果代表的是站點所在的斷面水體的垂向平均鹽度超標小時數(shù)除以1 a小時數(shù)(365 d×24 h)的算術(shù)結(jié)果，實測值則是某個取水口或者水閘所在位置的某一水深處的檢測值的統(tǒng)計結(jié)果，而在磨刀門水道的同一個斷面的不同位置(如主槽和邊灘、支汊)，或者在同一個位置的不同水深處，其鹽度值都有著較大的變化，因此，模型計算結(jié)果和實測統(tǒng)計值之間會存在一定范圍內(nèi)的誤差；③ 上游流量邊界不完全一致。

圖6為磨刀門水道咸淡水比例實測值及計算值的沿程變化示意。由磨刀門水道咸淡水比例沿程變化可知：計算結(jié)果和實測資料統(tǒng)計結(jié)果均表現(xiàn)為從下游向上游咸淡水比例逐漸降低；而從不同水文條件下的咸淡水比例變化可知：計算結(jié)果和實測資料統(tǒng)計結(jié)果均表明咸淡水比例變化與上游來流量有一定的負相關(guān)關(guān)系，枯水平均流量保證率越高則相應(yīng)的年咸淡水比例也越大。

圖6 磨刀門水道咸淡水比例實測值、計算值沿程變化示意

3.1.2東四口門

表3為東四口門黃埔站、中大站、三沙口站、南沙站和馮馬廟站的年淡水比例計算結(jié)果，圖7為珠江獅子洋水道咸水比例實測值與計算值的沿程變化。根據(jù)與實測資料的統(tǒng)計結(jié)果對比，三沙口站、南沙站和馮馬廟站的咸水比例計算與實測值之間的差異較小，如三沙口站的咸水比例模型計算結(jié)果為26%～37%，實測統(tǒng)計結(jié)果為33%～38%；南沙站的模型計算結(jié)果為33%～45%，實測統(tǒng)計結(jié)果為35%～40%；馮馬廟站淡水比例計算結(jié)果為12%～18%，實測統(tǒng)計結(jié)果為25%。黃埔站和中大站的咸水比例模型計算結(jié)果要小于統(tǒng)計結(jié)果，而且誤差較大，與實測資料的局限和統(tǒng)計方法有關(guān)。該站實測資料的咸淡水比例計算方法是采用每日高低潮特征時刻鹽度值計算，1 a中1 d高低潮2次鹽度監(jiān)測資料，超標1次為0.5 d、超標2次為 1 d，咸水比例=含氯度大于250 mg/L /365 d。由此可知：三沙口站、南沙站和馮馬廟站的咸水比例計算結(jié)果與實測資料統(tǒng)計結(jié)果相符，而黃埔站和中大站的咸水比例模型計算結(jié)果與實測資料統(tǒng)計值存在較大誤差是合理的。因為三沙口站、南沙站和馮馬廟站位于河口咸水區(qū)的中段位置，在枯水期的大部分時間內(nèi)都不在咸水界移動范圍內(nèi)，因此，按超標1次為0.5 d、超標2次為1 d是比較符合實際的；而黃埔站及中大站位于河口咸水區(qū)的上段位置，枯水期的大部分時間內(nèi)都在咸水界移動范圍內(nèi)，假如其1 d鹽度超標1次，極有可能該天咸水界剛好在此位置，而導致實際其也只超標1～2 h，卻按0.5 d(12 h)計，最終咸水比例的統(tǒng)計結(jié)果將會比實際大很多。由此可知：模型計算得到的年咸水比例比較符合實際。珠江咸水比例從沿程變化可知：計算結(jié)果和實測資料統(tǒng)計結(jié)果均表現(xiàn)為從下游向上游咸水比例逐漸降低。

表3 珠江監(jiān)測站點咸淡水比例計算結(jié)果與實測結(jié)果 %

圖7 珠江咸淡水比例實測值、計算值沿程變化示意

3.2 典型年水文條件下珠江河口咸淡水比例的空間分布特征

圖8為計算得到的3個典型水文年條件下的珠江河口區(qū)全年咸淡水比例等值線分布示意。從圖8中可知：不同保證率典型年水文條件下的咸淡水比例等值線的分布規(guī)律基本一致。咸水比例0%等值線系計算典型年中咸水上溯的最遠距離，其位置主要受年最枯流量的影響，3組典型年水文條件下，咸水比例0%等值線的位置變化不大。從咸水比例等值線的疏密可知：0%～10%等值線較稀疏，且咸淡水比例越小各等值線之間的距離就越大，主要是咸水最大上溯距離在枯水期對流量的變化最為敏感(在較枯流量范圍內(nèi)，流量的較小變化能引起咸水界的較大距離的變化)；同時也因最枯流量等級流量出現(xiàn)的頻率很小。伶仃洋的咸水比例等值線分布與其鹽度場的分布一樣具有“東高西低”的特征，咸水比例100%的等值線在淇澳島與內(nèi)伶仃島之間。磨刀門河口咸水比例100%的等值線位于交杯灘外2～3 km處，形狀由河口內(nèi)向外突出。黃茅海的咸水比例100%等值線向西南上游凸起，咸水比例呈“西高東低”之勢，主要是由于徑流經(jīng)東部深槽下泄，減小了東部的咸水比例。

對比不同典型水文年下的咸淡水比例情況，可知不同典型年水文條件下相同的咸淡水比例等值線的位置變化不會太大，且一般情況下均是枯季平均流量保證率越高，各咸淡水比例等值線越往上游。不同頻率水文年下的咸淡水比例等值線的位置變化不足以說明珠江網(wǎng)河區(qū)的咸淡水分布具有相對穩(wěn)定性，在非極端水文條件下，網(wǎng)河區(qū)各處咸淡水比例均在一定的幅度內(nèi)變化。從圖8中亦可知：網(wǎng)河區(qū)偶有局部地區(qū)在枯季平均流量保證率高時候咸淡水比例等值線卻更往下游，主要是因為典型年的年內(nèi)流量分布并不一致，同時也從另一方面反映了珠江網(wǎng)河區(qū)的徑潮動力的復雜性。圖9為不同頻率水文年高要+石角水文站年內(nèi)日均流量分布，從圖9中知：枯季平均流量保證率90%的水文年的年內(nèi)日均流量并不都小于枯季平均流量保證率75%的水文年的年內(nèi)日均流量，如從0到100 d，保證率90%的日均流量便大于保證率為75%的日均流量，也從而導致網(wǎng)河區(qū)咸淡水比例等值線位置分布與對應(yīng)枯季平均流量保證率不完全一致。

圖8 枯季平均流量不同保證率水文條件下的咸淡水比例示意

圖9 不同頻率水文年高+石日均流量年內(nèi)分布示意

3.3 河口區(qū)主要取水口咸淡水比例

統(tǒng)計珠江河口區(qū)23家主要取水戶的不同保證率下的咸水比例可知，基本上所有取水戶的咸水比例均是在90%保證率的水文條件下最大，咸水比例受流量大小影響比較明顯，越是枯水年，各取水口的咸水比例越大。從河道沿程的取水口的年咸水比例變化來看，越往下游其咸水比例越大，而且越往下游，咸水比例增大的速度越快。其中，枯水期平均流量為50%、75%和90%保證率的水文條件下，取水口年咸水比例超過30%的取水戶分別為8戶、10戶和11戶。

將2019年珠江河口區(qū)主要電廠實際取水量及咸淡水比例折算應(yīng)計入用水總量控制部分水量，珠江河口主要電廠2019年實際取水量為319 807萬m3，按50%、75%和90%來水保證率計，取用咸淡水量分別為76 829萬m3、97 920萬m3和109 816萬m3。

4 結(jié)語

1) 建立基于非結(jié)構(gòu)網(wǎng)格的珠江河口三維鹽水數(shù)學模型，模型主要選取了2001年2月的珠江三角洲河網(wǎng)區(qū)同步測量水文資料進行潮位、流速、鹽度的率定。選擇2005年1月、2008年1月、2009年12月和2011年1月共4組鹽度實測資料對模型的鹽度計算進行驗證。水動力及鹽度驗證結(jié)果良好，能夠較好地反映原體實際情況。

2) 利用珠江河口三維鹽水數(shù)學模型，選取上游50%、75%、90%枯季來水保證率典型年，下游選取河口典型潮位過程，計算河口區(qū)不同水域(取水口)全年鹽度超標小時數(shù)并統(tǒng)計作圖，實現(xiàn)不同水文條件下珠江河口各水域鹽度值及不同水域鹽度超標時間的快速簡易查詢。提取珠江河口主要取水口咸淡水比例的計算結(jié)果與實測數(shù)據(jù)進行對比，發(fā)現(xiàn)計算結(jié)果與實測值基本相符。

3) 根據(jù)珠江河口主要取水戶咸淡水比例結(jié)果，在較典型水文條件下，取水口年咸水比例超過30%的取水戶超過8戶；主要電廠實際取水量中，咸淡水比例約為24.0%～34.3%。

4) 河口咸淡水區(qū)水資源利用具有獨特的區(qū)域特點，因此，建議利用本文提出的咸淡水比例計算方法，對河口大型取水戶的取水中的淡水和咸淡水進行量化，并將其量化結(jié)果應(yīng)用于河口區(qū)用水總量控制、水資源管理日常統(tǒng)計、水資源公報統(tǒng)計、水資源征收等方面，從而推進珠江河口區(qū)水資源精細化管理，為最嚴格水資源管理制度提供科學支撐。