磨刀門鹽水入侵特性及其動力機制研究

呂紫君，馮佳佳，孔 俊，羅照陽，潘明婕

(1.河海大學 江蘇省海岸海洋資源開發與環境安全重點試驗室，南京 210098;2.南京市水利規劃設計院股份有限公司，南京 210006)

鹽水入侵是沿海河口地區特有的一種水文現象，多發于枯水期，是河口研究中關鍵問題之一。國外關于鹽水入侵的研究從20世紀50年代開始，早期研究工作主要集中在河口環流、鹽水入侵長度、鹽淡水混合類型、最大渾濁帶等[1-3]。河口區的鹽淡水混合與輸運是一個極其復雜的過程，鹽水入侵受向陸和向海兩種鹽度輸送驅動力共同控制；Hansen[2]較早開始采用機制分解方法研究河口鹽度輸運，提出以河口流速及鹽度的垂、橫向變化來描述河口的鹽度通量；Lerczak[4]根據實測資料研究了哈德遜河口的鹽度輸運，將縱向鹽度輸運分成平流輸運、穩定剪切輸運及潮汐震動輸運。國內關于鹽水入侵的研究從20世紀80年代開始，主要集中在長江口[5]和珠江口[6]。

珠江口磨刀門水道沿途分布著眾多取水口，是珠海、中山、澳門等城市的主要飲用水水源地。21世紀以來，受航道疏浚、口門挖沙以及上游徑流量變化的影響，磨刀門鹽水入侵日趨嚴重，眾多學者針對這一問題開展了大量研究。賈良文等[7]對磨刀門水流及鹽度的變化過程進行分析，得出磨刀門枯季鹽淡水類型屬于緩混合型；韓志遠等[8]基于不同年份的潮位和地形資料，分析了磨刀門鹽水入侵加劇的原因；盧陳等[9]采用物理模型實驗方法對不同潮差驅動下鹽水入侵距離進行研究，結果表明存在潮差臨界值使得鹽水入侵距離最短。由于實地測量的時空覆蓋有限，而物理試驗的研究成本高，所以數值模擬已廣泛應用于河口鹽水入侵的研究，Gong和Shen[10]采用ELCIRC模型對磨刀門鹽水入侵與徑流和潮差的響應法則進行研究，表明徑流量的大小是影響鹽水入侵的關鍵因素；Wang等[11]、陳文龍等[12]利用FVCOM模型對磨刀門水道的鹽水入侵規律進行了分析，發現鹽水入侵最遠出現在小潮之后的中潮。

上述成果多為對磨刀門鹽水入侵規律與影響因素的研究，相應的物質輸運通量及其分解的研究偏少，鹽水入侵半月周期變化的動力學機制尚不明確，同時磨刀門鹽度數值模擬的精度及效率有待進一步提高。鑒于此，本文基于一種最新發展的河口海洋模型SCHISM[13]，構建磨刀門水道三維水流鹽度數值模型，利用該模型對磨刀門鹽水入侵特性、縱向鹽度輸運進行分析，探究磨刀門鹽水入侵的動力機制。

1 數值模型建立

1.1 數值模型

研究中采用的模型為SCHISM[13]，該模型是基于SELFE[14]開發出的一個跨尺度湖泊-河流-河口-海洋數值模型，計算區域水平方向上采用三角形和四邊形混合網格，可以靈活地擬合復雜的地形和岸線；垂向上采用傳統的混合SZ坐標或者新增的LSC2坐標[15]，混合SZ坐標不僅能適應復雜不規則地形，而且能有效改善模型垂向的空間分辨率，克服河口區因地形變化大而引起的模擬水平梯度失真的難題；而LSC2坐標在保證深海區垂向網格精度的同時，避免淺水區垂向層數過多。模型采用高階ELM方法處理動量方程中的對流項，采用迎風格式或隱式二階TVD格式處理物質輸運方程中的對流項，既可提高計算的時間步長，又可保證模型計算穩定性和效率。采用GLS模型求解摻混系數，并提供了與通用的GOTM模塊的接口。關于SCHISM模型的其他詳細介紹，可見參考文獻[13]，模型在水動力模型基礎上耦合了波浪、生態、水質、風暴潮等模塊，已經用于全球的水動力相關問題的研究與預報，如北海和波羅海[16]等。

1.2 計算區域與網格

磨刀門三維數值模型共有2個開邊界，上游開邊界位于百頃頭以上的外海大橋附近，給流量過程，下游至磨刀門口外-15 m等深線附近，給水位過程；邊界流量、水位及鹽度均從已建的珠江口二維模型中提取，穩定的初始鹽度場通過循環運行模型得出。模型水平方向采用非結構三角形網格，網格共23 283個，節點共13 168個；模型垂向采用混合SZ坐標，平均海平面45 m以下采用Z坐標，分1層；45 m以上采用S坐標，分11層，每層厚度隨地形變化。網格分辨率的大小直接影響著模型的計算效率和精度，故分辨率分區域給定，口門區域分辨率取為1 500 m，而為確保鹽度模擬精度，河道內分辨率取在40 m以下。

1.3 模型計算參數

模型計算時間為2009年1月2日～1月31日，共30 d，時間步長為150 s；模型全域給恒定曼寧系數，大小取0.015；模型為斜壓模型，采用冷啟動，不考慮底床變形和風應力，采用GLS湍流混合模型求解摻混系數，采用隱式TVD2格式處理輸運方程中的對流項，采用高階ELM方法處理動量方程中的對流項，干濕網格判斷最小水深均為0.1 m。

圖1 磨刀門水深地形及測點位置圖Fig.1 Topography of the Modaomen estuary and the locations of measurement stations

2 模型驗證

采用2009年1月的實測資料對建立的數值模型進行驗證，驗證點位置分布見圖1，其中三灶站、燈籠山站、竹銀站為水位測站，M1 ～ M8點為流速和鹽度測點。為定量評估驗證精度，引入均方根誤差(RMSE)，其表達方式如下

(1)

式中：Xmod為模型計算值；Xobs為實際測量值；N為觀測次數。鑒于篇幅所限，本文只列出了部分驗證成果，如圖2～圖4所示，圓點為實測值，實線為模擬值。驗證表明，計算值與實測值吻合較好，誤差均在允許范圍之內。總體上，模型能較好地反演磨刀門水流及鹽度的變化過程，因此建立的三維水流鹽度數值模型可用于鹽水入侵機制的研究。

圖3 M2、M6測點流速驗證Fig.3 Calibration of flow velocity at M2 and M6 station

圖4 M2、M6測點鹽度驗證Fig.4 Calibration of salinity at M2 and M6 stations

3 結果分析

3.1 磨刀門鹽水入侵特征

圖5 M6點表、底層鹽度時間過程線Fig.5 Temporal variation of surface and bottom salinity at M6 station

珠江河口的潮汐類型為不規則半日潮，一天中出現2次高潮和2次低潮，潮差和潮歷時明顯不等。圖5給出了M6點表、底層的鹽度時間變化過程，為說明鹽度隨大小潮潮型的變化，同時給出了三灶站的潮位過程。由圖可知，M6點的底層鹽度變化在0～9 ppt之間，表層鹽度變化在0～6 ppt之間；且鹽水入侵過程有著明顯的半月周期變化規律，表底層鹽度在小潮后期驟然升高，并在中潮時達到峰值，之后隨著潮差增大，鹽度慢慢變小。

圖6和圖7分別給出了小潮、大潮期間磨刀門水道潮平均的表、底層的流場和鹽度場分布情況，小潮期間(圖6)，表層余流在整個磨刀門水道包括洪灣水道及鶴州水道都是向海運動，磨刀門水道流速有明顯的橫向差異，東側流速明顯較西側大，這是因為東側為深槽，是重要潮流通道，流速相對較大，而西側為淺灘，底摩阻相對較大，削弱了潮流動力，流速相對較??；底層余流在磨刀門水道主要向陸運動，在西岸淺灘區域則向海運動，在洪灣水道上段向海運動，下段則向陸運動；鹽水沿磨刀門水道東側深槽入侵，底部10 ppt鹽水可入侵到掛定角以上4 km附近，表層10 ppt等鹽度線則處于口門大橫琴附近，表底層鹽度差異大，鹽水分層明顯。大潮期間(圖7)，表層余流在磨刀門口外海區、口內河道和洪灣水道仍是向海運動，底層余流在磨刀門水道深槽區向陸運動，在淺灘則向海運動；與小潮相比，磨刀門水道內鹽水入侵減弱，底部10 ppt鹽水回落到掛定角以上1.5 km處，鹽水分層較弱。

從潮周期平均的流場和鹽度場分布可以看出，磨刀門水道東側深槽是高鹽水入侵磨刀門的重要通道，小潮期間磨刀門水道攔門沙至掛定角段底部積聚較高的鹽水，隨著潮差增大，積聚的鹽水被輸送至上游。

圖6 小潮期間表層(左)及底層(右)流場、鹽度分布Fig.6Surface(left)andbottom(right)meancurrentandtidally-averagedsalinityatneaptide圖7 大潮期間表層(左)及底層(右)流場、鹽度分布Fig.7Surface(left)andbottom(right)meancurrentandtidally-averagedsalinityatspringtide

3.2 磨刀門鹽水入侵動力機制分析

3.2.1 鹽淡水混合特征

為了進一步分析磨刀門鹽水入侵的動力機制，利用數值模型計算結果，繪制了漲急、落急時刻A-A縱斷面的瞬時鹽度場，縱斷面位置如圖1所示。選取分層系數(n)來描述不同時刻鹽度的分層特征，n為測點表、底層的鹽度差與測點垂線平均鹽度的比值，定義為

(2)

圖8-a為小潮期間漲急、落急時刻磨刀門縱斷面鹽度分布圖，圖上橫軸為與外海A點的距離，流速箭頭代表沿河道的流速(下同)。小潮期間，三灶站的漲潮潮差為0.93 m，落潮潮差為0.68 m，潮汐動力較小，縱斷面的鹽度等值線走勢比較平緩，沿程鹽度分層明顯。漲急時刻，口外底部水體鹽度變大，最高鹽度可達30 ppt，而表層鹽度值較低，隨著漲潮動力增強，口外高鹽水團沿河道底部向上游推進，表層水體受其擠壓作用，形成明顯的高、低鹽水分界面；此時，n2= 0.77，n4= 1.65，n6= 0.01，磨刀門水道的鹽水幾乎處于高度分層狀態，且越往上游，分層越明顯，直到上游燈籠山附近，鹽水濃度很低，鹽度均勻分布，分層現象消失。落急時刻，中層及表層鹽水隨落潮流退至外海，但口門至洪灣水道段底部仍積蓄著高鹽水；此時，n2= 0.47，n4= 0.91，n6= 1.28，水道下游的鹽水處于緩混合狀態，而洪灣水道至燈籠山的鹽水仍處于高度分層狀態。

中潮期間(圖8-b)，三灶站的漲潮潮差為1.01 m，落潮潮差為1.73 m，與小潮期間相比，縱斷面的鹽度等值線趨勢有所傾斜，鹽水分層減弱。漲急時刻，口外鹽水楔進入河道，與表層低鹽水混合后整體向上游入侵，水道內鹽水濃度明顯變高，鹽水入侵距離增加；此時，n2= 0.51，n4= 0.98，n6= 1.46，在鹽水楔鋒面所在范圍內，鹽水處于高度分層狀態。落急時刻，口外高鹽水稍有退卻，鹽水入侵距離達到峰值；此時，n2= 0.34，n4= 0.55，n6= 0.35，鹽淡水混合充分，河道內的鹽水基本處于緩混合狀態。

大潮期間(圖8-c)，三灶站的漲潮潮差為1.31 m，落潮潮差為2.01 m。漲急時刻，口外高鹽水團進入河道，相比中潮漲急時刻，磨刀門水道內鹽水濃度變低，鹽水入侵距離減??；此時，n2= 0.62，n4= 1.24，n6= 0.70，磨刀門水道的鹽水分層加強。落急時刻，口外高鹽水隨落潮流退出河道，此時，n2= 0.41，n4= 0.89，n6= 0.67，鹽水基本處于緩混合狀態，大潮潮汐動力強，一天中鹽水濃度變化劇烈，鹽水“大進大出”。

3.2.2 鹽通量分解

本文采用Lerczak等[4]提出的計算鹽度通量的方法，將河口的鹽度輸運分為平流輸運(advection transport)、穩定剪切輸運(steady shear transport)、潮汐震蕩輸運(tidal oscillatory transport)三部分。斷面總鹽通量可以表示如下

FS==<(u0+uE+uT)(S0+SE+ST)dA>

(3)

式中：u為軸向流速，m/s；s為鹽度，ppt；A為橫斷面面積，m2；下標0為潮平均的余流；下標E為剪切流；下標T為潮流波動。MacCready[17]認為公式中的交叉項可以近似為零，得出鹽度通量的分解表達式如下

FS=<(u0s0uESE+uTST)dA>=Qfs0+FE+FT

(4)

式中：FS為總鹽通量；Qfs0為平流輸運項；FE為穩定剪切項；FT為潮汐震蕩項。平流輸運項是由上游流量以及外海stokes漂流引起的，主要是將上游鹽水輸運至外海；穩定剪切項主要由垂向剪切流引起，而潮汐震蕩項是由流速與鹽度的潮內變化相互作用產生，它們則主要是將外海的高鹽水輸運至河口上游。

以口門處的M2點為起始點，向上游河道每隔5 km取一個橫斷面，共選取6個橫斷面，計算并分析各斷面的鹽通量在一個潮周期過程中的變化情況。圖9給出了斷面2及斷面5的鹽通量分量Qfs0、FE和FT以及總鹽通量FS在一個完整潮周期過程中的變化情況，圖中正值表示向海輸運；負值表示向陸輸運，斷面具體位置如圖8所示。

圖9-c為斷面2的鹽通量分量，在斷面2處，Qfs0始終為正值，即平流輸送作用向海輸送鹽水，且在小潮期間平流輸送作用弱，在大潮期間平流輸送作用強；FE始終為負值，即穩定剪切作用向陸輸送鹽水，且在小潮期間穩定剪切作用強，在大潮期間穩定剪切作用明顯減弱；潮汐震蕩作用FT在斷面2處向陸輸送鹽水，且在大潮期間最強，在小潮期間最弱；總鹽通量FS在小潮期間通過斷面2向河道內輸送鹽水，而在大潮期間向外海輸送鹽水。

圖9-d為斷面5的鹽通量分量，在斷面5處，穩定剪切輸送、平流輸送和潮汐震蕩的規律與斷面2基本相同，但輸送強度變小；FS在后半月周期的小潮期間向陸輸送鹽水，隨后在中潮轉大潮、大潮期間向海輸送鹽水，而在前半月周期(1月6日～1月17日)，FS基本為0，表明在Qfs0、FE和FT三者作用下，此期間磨刀門水道內的日潮平均總鹽量處于平衡狀態。

9-a 三灶站水位 9-b 鹽水入侵長度

9-c 斷面2的鹽通量分量 9-d 斷面5的鹽通量分量圖9 鹽通量分量過程圖Fig.9 Decomposition of salt transport flux

4 結論

基于SCHISM模型，對磨刀門水道鹽水入侵特性、縱向鹽度輸運進行了分析，得出以下結論：

(1)枯季磨刀門鹽水入侵一般特性為：小潮期間，磨刀門水道鹽水分層明顯，攔門沙至掛定角段底部積聚高鹽水，鹽水入侵距離較小；中潮期間，鹽水混合充分，小潮時積聚的底部高鹽水逐漸摻混至表層，并被輸送至上游，鹽水入侵距離最遠；大潮期間，一天中鹽水濃度變化劇烈，鹽水“大進大出”，鹽水入侵距離比中潮期間有所減小，但大于小潮期間的入侵距離。

(2)磨刀門日平均鹽度受平流輸送和穩定剪切輸送共同作用，鹽量輸運在大部分時間內處于不平衡狀態。小潮期間，總鹽通量向陸，穩定剪切作用強，平流輸運作用弱，鹽量在河道里累積；中潮期間，總鹽通量由向陸轉為向海，河道內的日平均鹽度最大，鹽水入侵距離最遠；大潮期間，總鹽度通量向海，驅動鹽水向陸輸送的穩定剪切和潮汐震蕩作用小于驅動鹽水向海輸送的平流輸送作用，鹽量被沖出河口。

參考文獻：

[1]PRITCHARD D W. Salinity distribution and circulation in the Chesapeake Bay estuarine system [J]. Journal of Marine Research, 1952, 11(2): 106-123.

[2]HANSEN D V, RATTRAY M. Gravitational circulation in straits and estuaries [J].J.mar.res,1965, 11(3): 319-326.

[3]SIMMONS H B, BROWN F R. Salinity effects on estuarine hydraulics and sedimentation [C]//Tanaka H. Process of 13th IAHR (3).Kyoto, Japan: International Association for HydroEnvironment Engineering and Research, 1969:192-216.

[4]LERCZAK J A, GEYER W R, CHANT R J. Mechanisms Driving the Time-Dependent Salt Flux in a Partially Stratified Estuary [J]. Journal of Physical Oceanography, 2006, 36(12): 2 296-2 311.

[5]CHEN W, CHEN K, KUANG C, et al. Influence of sea level rise on saline water intrusion in the Yangtze River Estuary, China [J]. Applied Ocean Research, 2016, 54:12-25.

[6]歐素英. 珠江三角洲咸潮活動的空間差異性分析 [J]. 地理科學, 2009, 29(1): 89-92.

OU S Y. Spatial difference about activity of saline water intrusion in the Pearl River Delta [J]. Scientia Geographica Sinica, 2009,29(1): 89-92.

[7]賈良文, 吳超羽, 任杰, 等. 珠江口磨刀門枯季水文特征及河口動力過程 [J]. 水科學進展, 2006, 17(1): 82-88.

JIA L W, WU C Y, REN J, et al. Hydrologic characteristics and estuarine dynamic process during the dry season in Modaomen Estuary of the Pearl River [J]. Advances in water science, 2006, 17(1): 82-88.

[8]韓志遠, 田向平, 劉峰. 珠江磨刀門水道咸潮上溯加劇的原因 [J]. 海洋學研究, 2010, 28(2): 52-59.

HAN Z Y, TIAN X P, LIU F. Study on the causes of intensified saline water intrusion into Modaomen estuary of Pearl River in recent years [J]. Journal of Marine Science, 2010, 28(2): 52-59.

[9]GONG W, SHEN J. The response of salt intrusion to changes in river discharge and tidal mixing during the dry season in the Modaomen Estuary, China [J]. Continental Shelf Research, 2011, 31(7): 769-788.

[10]盧陳, 袁麗蓉, 高時友, 等. 潮汐強度與咸潮上溯距離試驗 [J]. 水科學進展, 2013, 24(2): 251-257.

LU C,YUAN L R,GAO S Y, et al. Experimental study on the relationship between tide strength and salt intrusion length [J]. Advances in water science, 2013, 24(2): 251-257.

[11]WANG B, ZHU J, WU H, et al. Dynamics of saltwater intrusion in the Modaomen Waterway of the Pearl River Estuary [J]. Science China Earth Sciences, 2012, 55(11): 1 901-1 918.

[12]陳文龍, 鄒華志, 董延軍. 磨刀門水道咸潮上溯動力特性分析 [J]. 水科學進展, 2014, 25(5): 713-723.

CHEN W L, ZOU H Z, DONG Y J. Hydrodynamic of saltwater intrusion in the Modaomen waterway [J]. Advances in water science, 2014, 25(5): 713-723.

[13]ZHANG Y J, YE F, STANEV E V, et al. Seamless cross-scale modeling with SCHISM [J]. Ocean Modelling, 2016, 102:64-81.

[14]ZHANG Y, BAPTISTA A M. SELFE: a semi-implicit Eulerian-Lagrangian finite-element model for cross-scale ocean circulation [J]. Ocean modelling, 2008, 21(3): 71-96.

[15]ZHANG Y J, ATELJEVICH E, YU H C, et al. A new vertical coordinate system for a 3D unstructured-grid model [J]. Ocean Modelling, 2015, 85:16-31.

[16]ZHANG Y J, STANEV E V, GRASHORN S. Unstructured-grid model for the North Sea and Baltic Sea: validation against observations [J]. Ocean Modelling, 2016, 97: 91-108.

[17]MACCREADY P, GEYER W R. Estuarine salt flux through an isohaline surface [J]. Journal of Geophysical Research Atmospheres, 2001, 106(C6): 11 629-11 638.