林若蘭，卓文珊，曾 珂，梁廖逢，劉 茹，劉祖發(fā)

(1.中山大學(xué)地理科學(xué)與規(guī)劃學(xué)院，廣東 廣州 510275； 2.中山大學(xué)水資源與環(huán)境研究中心，廣東 廣州 510275；3.中山大學(xué)測試中心, 廣東 廣州 510275； 4.華南地區(qū)水循環(huán)與水安全廣東省普通高校重點實驗室，廣東 廣州 510275)

河口咸潮入侵對珠三角城市群生產(chǎn)、生活、生態(tài)各方面用水有著重要的影響[1]。目前針對珠江口的咸潮問題已經(jīng)進行了較多的研究，從方法上主要分為3類：基于實測資料的數(shù)學(xué)統(tǒng)計分析[2-7]、基于數(shù)值模型的分析與計算[8-11]及物理模型試驗[12]。但這些研究主要局限在徑流、潮流、海平面上升等方面對咸潮的影響進行研究，對風(fēng)這一影響因素的研究較少。實際上河口地區(qū)普遍存在著表層向海、底層向陸的環(huán)流，除了密度梯度效應(yīng)外，風(fēng)是形成這種環(huán)流的另一重要驅(qū)動力。Hansen等[13]研究表明風(fēng)由河道上游向海吹時加強河口兩層環(huán)流，反方向則抑制環(huán)流發(fā)展；Weisberg等[14-16]也發(fā)現(xiàn)了風(fēng)作用下的兩層流結(jié)構(gòu)；Manca等[17]對臺伯河河口的研究發(fā)現(xiàn)風(fēng)速大于4 m/s時能控制表層海流并抬升潮位，從而解釋臺伯河河口在徑流量較大的情況下鹽水入侵仍能發(fā)生；Héctor等[18]通過實測水文資料對弱潮汐型河口Jamapa進行研究，發(fā)現(xiàn)下河口風(fēng)(局部作用)和西北風(fēng)(遠距離作用)都能增強鹽水楔入侵，影響鹽水楔位置；秦志新等[19]研究發(fā)現(xiàn)在西北偏北的大風(fēng)作用下，沐官島海灣水庫突然出現(xiàn)垂向混合從而導(dǎo)致表層鹽度超標的風(fēng)險較大。由此可見，風(fēng)對鹽水入侵的影響不容忽視。針對珠江口也有學(xué)者進行了相關(guān)研究，包蕓等[20]基于Backhaus的三維斜壓非線性模型論述了均勻西南風(fēng)場對斜壓場的影響；劉歡等[21]運用Ecomsed三維斜壓模型計算風(fēng)對垂向環(huán)流形成的影響，推斷枯季期間磨刀門水道在東北季風(fēng)影響下形成環(huán)流的情況。但對不同風(fēng)向下的潮流運動和鹽度分布進行對比的研究較少。李越等[22]對枯季不同風(fēng)向?qū)V州附近水道咸潮的影響進行分析，但采用一維、二維聯(lián)解數(shù)學(xué)模型，無法反映表底層水體對風(fēng)的響應(yīng)差異。

本文根據(jù)珠江東四口門枯季實測資料，采用MIKE3建立三維水動力模型，在驗證模型的可靠性基礎(chǔ)上，分析在無風(fēng)、北風(fēng)、東北風(fēng)和東風(fēng)的作用下枯水期東四口門的潮流運動和鹽度分布的變化，對咸潮入侵問題的認識起補充作用，旨在為河口生態(tài)環(huán)境的維護、水資源管理等提供參考。

1 模型的建立

1.1 控制方程及計算方法

MIKE3數(shù)值模型可模擬具有自由表面的三維流動系統(tǒng)，同時考慮質(zhì)量守恒、動量守恒、鹽度及溫度守恒，局部鹽度、溫度、壓力平衡。其數(shù)學(xué)基礎(chǔ)是雷諾平均化的N-S方程，包括紊流影響、密度變化及溫鹽平衡方程。水動力學(xué)模塊采用交替方向隱式迭代法(ADI方法)對質(zhì)量及動量守恒方程進行積分，且對產(chǎn)生的數(shù)學(xué)矩陣采用雙精度掃描法進行求解[23]。

1.2 計算區(qū)域及網(wǎng)格劃分

計算區(qū)域為北江三角洲水系、東江三角洲水系及伶仃洋河口灣，上游到達博羅站、老鴉崗站、三水站以及南華站，下游到達伶仃洋25 m等深線附近(圖1)。東四口門中虎門對應(yīng)水文站點為大虎站、泗盛圍站，蕉門、洪奇門、橫門對應(yīng)水文站點分別為南沙站、馮馬廟站、橫門站。上游水動力邊界條件為水位或者流量，下邊界條件采用經(jīng)過校正的全球潮汐預(yù)測模型生成同步的外海邊界的潮位。邊界鹽度條件采用實測鹽度值，其中三水、南華、博羅實測鹽度為0，老鴉崗的邊界鹽度由非逐時數(shù)據(jù)插值而成，外海開邊界鹽度值設(shè)為35‰。初始鹽度場由研究區(qū)18個實測站點鹽度數(shù)據(jù)插值而成，有利于計算穩(wěn)定并加快收斂，率定驗證時的風(fēng)條件采用同步的深圳寶安機場的數(shù)據(jù)。為更好地擬合多變的河口灣地形和曲折岸線，采用非結(jié)構(gòu)化三角網(wǎng)格，外海開邊界處采用5 km的分辨率，內(nèi)河道采用500 m左右的分辨率，部分細河道采用50～100 m左右的分辨率，計算網(wǎng)格共包含33 445個網(wǎng)格，21 766個節(jié)點，垂向上采用σ坐標等距劃分5層,采用1999年的河道地形數(shù)據(jù)。模型采用斜壓模式，考慮溫度和鹽度對密度的影響，密度計算根據(jù)UNESCO海水標準方程。設(shè)置干濕邊界的臨界判斷值，對河床糙率等重要參數(shù)進行率定，水平擴散性系數(shù)采用Smagorinsky公式進行計算。模型考慮徑流、科氏力、潮流、風(fēng)等重要動力因子的影響，不考慮降水、蒸發(fā)以及水工結(jié)構(gòu)物等因素的影響。

圖1 模型范圍內(nèi)河道地形及測站分布

2 模型的驗證

近幾十年來人類采砂、疏浚等活動對河口地形產(chǎn)生較大影響，為保證地形和水文資料相匹配，采用與地形年份最接近的珠江口2001年2月8—16日的水位、流量、鹽度數(shù)據(jù)進行率定和驗證。其中水位Z和流量Q為逐時數(shù)據(jù)，鹽度S為非逐時數(shù)據(jù)。僅列出東四口門附近站點的模擬情況，模擬結(jié)果與實測值的對比情況見圖2、圖3、圖4。

(a) 大虎

(b) 南沙

(c) 馮馬廟

(d) 橫門

(a) 大虎

(b) 南沙

(c) 馮馬廟

(d) 橫門

(a) 泗盛圍

(b) 南沙

(c) 馮馬廟

(d) 橫門

珠江口潮汐類型為不規(guī)則半日潮，一天內(nèi)分別有兩個不等的高潮和低潮，由圖2可知各站點從大潮到小潮的模擬水位與實測水位值具有很好的擬合度，無論是振幅還是相位均比較吻合，相關(guān)系數(shù)R均大于0.9，納西效率系數(shù)NS均大于0.8。由圖3可知，各站點流量大小及漲落潮流向的變化也擬合較好，R均大于0.9，NS均大于0.7，但4個站點流量振幅整體偏小，可能是地形數(shù)據(jù)分辨率不足所致。由圖4可知，從大潮到小潮，鹽度逐漸降低，橫門、馮馬廟由于口外淺灘淤積，鹽水入侵較弱，總體鹽度較低，南沙漲潮和落潮時鹽度相差較大，漲潮時鹽水從蕉門的龍穴南水道和鳧洲水道兩個通道進入蕉門，鹽度增加，落潮時較大的下泄徑流量使鹽度大大降低。由于虎門是潮汐型河口，鹽水上溯較強，泗盛圍的鹽度整體較大，即使小潮落潮時也大于2‰，模型基本能反映鹽度的變化規(guī)律。

3 模型的應(yīng)用

風(fēng)是影響潮汐運動的驅(qū)動力之一，在平面上，風(fēng)通過剪切應(yīng)力引起水體的平面輸送；垂向上，風(fēng)通過垂向氣壓的變化促使表底層水體交換進而影響鹽度的分布，影響河口環(huán)流。為了進一步探究風(fēng)的影響，應(yīng)用經(jīng)過率定、驗證的模型，通過數(shù)值試驗來分析風(fēng)對東四口門咸潮的影響。珠江口枯季常風(fēng)向為北，次常風(fēng)向為東、東北，2000年12月至次年2月季風(fēng)平均風(fēng)速為4.9 m/s[24]，因此選取北風(fēng)、東北風(fēng)和東風(fēng)3種風(fēng)向來進行研究，風(fēng)速統(tǒng)一為5 m/s，水文條件選取2001年2月至4月的日平均流量，對 2001年2月7日至4月1日進行模擬，模擬時段包含4個大潮和4個小潮，分析不同風(fēng)作用下流速、鹽度、潮通量的變化。風(fēng)對水體的表面摩擦力取決于風(fēng)速，剪力值τs計算公式為

τs=ρacd|μw|μw

(1)

式中：ρa為空氣密度；μw為海平面上10 m處所測量到的風(fēng)速；cd為拖曳系數(shù)，根據(jù)Wu[25]的經(jīng)驗公式進行計算。

4 結(jié)果分析

4.1 流速

虎門鹽水楔入侵距離較長。以虎門為例，分別計算模擬時段大小潮期大虎站各層水體在北風(fēng)、東北風(fēng)、東風(fēng)3種條件下的漲落潮流速及平均流速變化率(表1)。大潮時，潮汐動力較強，風(fēng)的影響有限，各層水體的漲落潮流速變幅在1%以內(nèi)，平均流速變幅在11%以內(nèi)，由于底層流速較小，所以底層流速變幅較大。落潮時，北風(fēng)通過表面應(yīng)力使表層落潮流速增大(相比于無風(fēng)條件)，底層產(chǎn)生向上游的補償流，使底層落潮流速減小。而東風(fēng)對落潮流速的影響則相反，抑制表層落潮流，促進底層落潮流。漲潮時，北風(fēng)使表層漲潮流速減小，底層漲潮流速增大，東風(fēng)則相反，促進表層漲潮流，抑制底層漲潮流。東北風(fēng)的影響與北風(fēng)相同，但由于風(fēng)向與虎門水道呈一定夾角，影響較小。大潮時虎門各層水體的平均流速向海，表層比底層大，在北風(fēng)作用下，中底層的向海平均流速減小，底層減小10.8%，而表層的向海平均流速增大5.3%。東北風(fēng)作用下底層平均流速減小2.7%，表層平均流速增大1.9%。東風(fēng)作用下表層的向海平均流速減小2.8%，底層平均流速增大7.0%。

表1 在不同風(fēng)向下大潮、小潮期間大虎站各層水體的流速變化率

注：1層為底層；2、3、4層為中層；5層為表層。

小潮時流速對風(fēng)向的響應(yīng)與大潮時相同(表1)，北風(fēng)使表層漲潮流速減小，落潮流速增大，并產(chǎn)生補償流，使中底層漲潮流速增大，落潮流速減小。東北風(fēng)作用下的流速變化與北風(fēng)相同但變化幅度較小，東風(fēng)則完全相反。由于小潮時潮汐動力減弱，風(fēng)力的影響相對增強，流速的變化幅度整體變大。小潮時在北風(fēng)作用下，中底層的向海平均流速減小，其中底層減小40.5%，表層的向海平均流速增大9.1%。東北風(fēng)作用下底層平均流速減小9.4%，表層平均流速增大2.6%。東風(fēng)作用下表層的向海平均流速減小6.4%，底層平均流速增大26.6%。

從垂向平均來看，北風(fēng)使垂向平均流速減小，東北風(fēng)也使其減小但幅度較小，東風(fēng)則能使其增大。此外，與大潮時相比，大虎站小潮時多次出現(xiàn)底層水體的周期平均流速向陸，這與小潮時鹽淡水混合作用弱，分層增強及密度流速度大有關(guān)[26]。用Hansen等[27]提出的無量綱環(huán)流系數(shù)G(潮周期平均的下泄流速與斷面平均流速的比值)來表征垂向環(huán)流強度，當G=1,表底層流向相同；當G>1，垂向環(huán)流形成。計算模擬時段小潮期間多個太陰日周期的G，在北風(fēng)作用下，小潮時G>1的頻次增加，且G增大；在東北風(fēng)的作用下，G>1的頻次不變，但G增大；在東風(fēng)作用下，G>1的頻次減少。結(jié)果表明，北風(fēng)、東北風(fēng)促進虎門垂向環(huán)流的發(fā)展，而東風(fēng)抑制垂向環(huán)流。

4.2 垂向鹽度

選取模擬時段第一個和最后一個(風(fēng)持續(xù)作用45 d后)大小潮期的鹽度進行分析,為方便描述用大潮a、小潮a、大潮b、小潮b表示。圖5為風(fēng)作用2 d后的大潮a時各風(fēng)況下大虎山南到舢板洲的周期平均鹽度縱剖面圖。在北風(fēng)作用下，虎門表層鹽度整體減小，底層鹽度增加，由于虎門潮汐動力較強，5 m/s的風(fēng)作用有限，鹽度等值線移動的距離不大。北風(fēng)產(chǎn)生的表面應(yīng)力的分力與口門落潮流方向相同，使表層落潮流增強，更多的表層淡水被帶到外海，表層鹽度略微下降，底層向上游的補償流增強，底層的鹽水楔入侵略微增強，底層鹽度比無風(fēng)時大。在東北風(fēng)作用下，鹽度分布與無風(fēng)時相近，這是由于東北風(fēng)方向與虎門水道有一定夾角，且虎門潮汐動力較強，對漲落潮影響有限。東風(fēng)使表層漲潮流速增大，落潮流速減小，促進表層鹽水向上游的輸運，表層鹽度比無風(fēng)時略大，由于中底層的向海流速增大，底層鹽水的楔入距離減小。從垂向平均鹽度來看，與無風(fēng)條件相比，北風(fēng)、東北風(fēng)使大虎山南到舢板洲的平均鹽度增大，且北風(fēng)的變幅更大，東風(fēng)使平均鹽度減小。

(a) 無風(fēng)

(b) 北風(fēng)

(c) 東北風(fēng)

(d) 東風(fēng)

從大潮a到小潮a(風(fēng)作用9 d后)，潮差減小，潮汐動力減弱，風(fēng)的作用對虎門鹽度分布的影響較明顯(圖6)。與無風(fēng)條件相比，小潮時，在北風(fēng)作用下鹽度垂向梯度增加，分層增強，虎門表底層鹽度增加，大于10‰的鹽水整體上溯1.6 km，楔入河口的底層鹽水鹽度達18‰。表層鹽度的增加是由于底層鹽水楔入侵增強較明顯，表層和底層發(fā)生垂向紊動混合，表層鹽度也略微增大。在東北風(fēng)的作用下，上游段表底層變化不大，但口門處底層16‰的鹽水楔入增強。在東風(fēng)作用下，垂向鹽度梯度明顯減小，分層減弱，楔入河口的鹽度明顯減小，降至14‰，垂向混合作用下表層鹽度也減小。小潮時垂向平均鹽度的變化規(guī)律與大潮時相同。從2月到3月底，隨著上游來水流量的增大，鹽水逐漸后退，同時隨著時間的積累，不同風(fēng)向下的鹽度分布差異較大(圖7和圖8)。北風(fēng)明顯抑制鹽水的后退，使鹽水在口門內(nèi)停留的時間較長，東北風(fēng)也有相似的作用但不明顯，東風(fēng)加快鹽水的后退。

(a) 無風(fēng)

(b) 北風(fēng)

(c) 東北風(fēng)

(d) 東風(fēng)

(a) 無風(fēng)

(b) 北風(fēng)

(c) 東北風(fēng)

(d) 東風(fēng)

4.3 潮通量

物質(zhì)的滯留、積蓄作用對河口海岸環(huán)境的影響至關(guān)重要，計算模擬時段典型大小潮期橫門、馮馬廟、南沙、大虎4個站點的潮通量變化(表2)，向陸為正，向海為負，表中的大潮、小潮的太陰日周期的凈通量都為負，均為凈泄量。對于漲潮量，北風(fēng)、東北風(fēng)作用下大潮和小潮時橫門、馮馬廟、南沙、大虎4個站點都減少，東風(fēng)作用下大潮和小潮時4個站點都增加，但對于落潮量各個站點的變化規(guī)律不同。橫門在北風(fēng)作用下落潮量增加，這與北風(fēng)促進口門外的落潮流有關(guān)，在東風(fēng)、東北風(fēng)作用下落潮量減少，且東風(fēng)下的降幅更大，這與橫門水道東西走向有關(guān)。洪奇門與橫門落潮量的變化規(guī)律相似，不同點在于洪奇門在東北風(fēng)作用下落潮量是增加的，這與水道走向不同有關(guān)。南沙站和大虎站在北風(fēng)、東北風(fēng)的作用下落潮量減少。由4.1可知風(fēng)對大虎站表底層的漲落潮流速的影響不同，表現(xiàn)為表層順風(fēng)流速增加，逆風(fēng)流速減小，底層逆風(fēng)流速增加，順風(fēng)流速減小。由于伶仃洋東部潮汐動力較強，鹽度較高的鹽水能到達大虎、南沙這兩個站點處，形成鹽淡水異重流，在北風(fēng)、東北風(fēng)作用下異重流增強，可能使落潮量反而較小，而東風(fēng)則使落潮量增加。而橫門站和馮馬廟站徑流作用遠大于潮汐動力，鹽度較高的鹽水團不能到達該位置處，平均鹽度小于1‰(實測數(shù)據(jù))，鹽淡水密度差小，不存在明顯異重流。由于大潮時潮汐動力較強，風(fēng)的影響相對較小，因此小潮時漲落潮量的變化幅度比較大，但也不超過5%。

(a) 無風(fēng)

(b) 北風(fēng)

(c) 東北風(fēng)

(d) 東風(fēng)

表2 各風(fēng)況下各站的潮通量變化

4.4 東四口門鹽水入侵

對風(fēng)作用45 d后的大小潮時期的表底層鹽度分布進行分析。圖9為無風(fēng)作用下大潮b時東四口門外的底層鹽度分布，可見鹽度分布存在橫向差異，東高西低，等鹽度線基本呈東北-西南走向。由于蕉門、洪奇門、橫門的徑流作用相對較強，虎門的徑流作用相對較弱，再加上科氏力使潮差向潮波傳播方向的右側(cè)遞增，等鹽度線向虎門凸出。淺灘和深槽的水動力特征存在差異，由于深槽區(qū)垂向環(huán)流發(fā)育，鹽水上溯更遠，等鹽度線在兩個深槽處凸出，呈現(xiàn)雙峰結(jié)構(gòu)，底層大于27‰的高鹽水上溯至內(nèi)伶仃島附近。在垂向上，由于海水密度較淡水密度大，海水從底部入侵，表層鹽度普遍比底層要低，由于表層水體受風(fēng)、波浪影響摻混作用強，表層水體等鹽度線的雙峰結(jié)構(gòu)不如底層明顯(圖10)?？菟谛〕逼陂g的鹽度分布與大潮相似，但由于小潮時潮汐動力較弱，鹽水上溯距離較小(圖略)。

由圖9可見，對于大潮時的底層水體鹽度分布，與無風(fēng)時相比，在北風(fēng)、東北風(fēng)的作用下，東四口門外底層水體鹽度整體增加，且增幅從東向西遞減，到伶仃洋西岸附近鹽度下降。在北風(fēng)的作用下，虎門口附近鹽度增幅最大，達4‰～6‰,這是由于虎門口鹽淡水交匯，表底層水體密度差大，且水深較大，垂向環(huán)流發(fā)育，而北風(fēng)促進垂向環(huán)流，使虎門口鹽水入侵加劇。到內(nèi)伶仃島以南徑流作用較小，表底層鹽度相差較小，分層減弱，北風(fēng)的影響減弱，鹽度增幅比虎門口門附近小。而伶仃洋西部為淺灘，在北風(fēng)摩擦應(yīng)力作用下，橫門、洪奇門向海潮通量增加(表2)，鹽度降低。在東北風(fēng)作用下，口門外鹽度變化整體與北風(fēng)相似，但和北風(fēng)相比虎門口附近鹽度增幅較小，伶仃洋西岸鹽度降幅較大。在東風(fēng)作用下，由于抑制河口垂向環(huán)流，虎門口門外底層鹽度下降，由于東風(fēng)使淡水流出蕉門、洪奇門、橫門后方向偏西南，因此越往西岸鹽度降幅越大，而礬石淺灘以南水域鹽度增加可能是由于東風(fēng)對鄰近海域的Ekman向岸運輸作用以及混合作用使鹽度略微增加。

(a) 無風(fēng)大潮b時的底層鹽度

(b) 北風(fēng)大潮b時的底層鹽度差

(c) 東北風(fēng)大潮b時的底層鹽度差

(d) 東風(fēng)大潮b時的底層鹽度差

(a) 無風(fēng)大潮b時的表層鹽度

(b) 北風(fēng)大潮b時的表層鹽度差

(c) 東北風(fēng)大潮b時的表層鹽度差

(d) 東北風(fēng)大潮b時的表層鹽度差

由圖10可見，對于大潮時的表層水體，由于在北風(fēng)、東北風(fēng)的作用下底層鹽水入侵加劇，表層水體和底層發(fā)生垂向混合，伶仃洋東部、中部表層水體的鹽度也增加。東風(fēng)的平流輸運作用促進中部表層水體與西部混合，使伶仃洋表層水體整體鹽度增加的面積大于東北風(fēng)、北風(fēng)。但對于虎門口表層水體，由于東風(fēng)抑制垂向環(huán)流，鹽水入侵較弱，口門內(nèi)表層、底層鹽度都下降。小潮時各風(fēng)作用下的表底層鹽度變化規(guī)律與大潮時相同。

結(jié)果表明，北風(fēng)、東北風(fēng)的持續(xù)作用使虎門鹽水入侵加劇，虎門口門處底層鹽度可增加4‰～6‰，表層鹽度增加2‰～4‰。由于鹽水可通過枕箱水道進入蕉門，對蕉門口鹽水入侵有一定影響。但北風(fēng)、東北風(fēng)對洪奇門、橫門的影響較小，表底層鹽度變化幅度不超過±2‰。東風(fēng)抑制東四口門鹽水入侵，東四口門外底層鹽度降低0～2‰，對于表層水體，虎門表層鹽度降低，其他三口門略微上升。

5 結(jié) 論

a. 表底層水體的流速對風(fēng)的響應(yīng)是相反的，表層水體順風(fēng)流動，底層水體逆風(fēng)流動，因此風(fēng)影響垂向環(huán)流的發(fā)展。在北風(fēng)作用下，虎門口小潮時垂向凈環(huán)流的出現(xiàn)頻次增加，且環(huán)流強度增大。在東北風(fēng)的作用下，垂向凈環(huán)流的頻次不變，但垂向環(huán)流強度增大。在東風(fēng)作用下，垂向凈環(huán)流的頻次減少。

b. 北風(fēng)、東北風(fēng)使東四口門漲潮量減少，東風(fēng)使?jié)q潮量增加。而落潮量的變化則不僅與風(fēng)向相關(guān)，還與河道走向、垂向環(huán)流有關(guān)，因此四個口門的落潮量對風(fēng)的響應(yīng)有所不同。

c. 不同風(fēng)向下河口灣鹽度的空間分布也有所不同，從而影響到上游口門內(nèi)的咸潮活動。風(fēng)的持續(xù)作用對虎門內(nèi)鹽水入侵影響較大，北風(fēng)、東北風(fēng)促進虎門鹽水的入侵，當上游來水流量增大時，北風(fēng)、東北風(fēng)抑制鹽水的后退，咸潮活動時間變長，東風(fēng)則相反。風(fēng)對于其他3個徑流型河口鹽水入侵有一定影響，但相對較小。因此，進行長時段的咸潮模擬預(yù)測時需要考慮季風(fēng)的影響，提高模擬預(yù)測精度。