強(qiáng)潮河口鹽水入侵特征數(shù)值模擬研究

路川藤，錢明霞*，夏威夷，丁偉，丁佩

（1.南京水利科學(xué)研究院，江蘇南京 210029；2.水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210098）

0 引言

入海河口在徑流和潮流的共同作用下，形成了特有的物理現(xiàn)象—鹽淡水混合。鹽淡水混合具有時(shí)空分布不均的特點(diǎn)，垂向上，水體可按鹽淡水的混合程度分為高度分層型、緩混合型和強(qiáng)混合型[1]；平面上，水體受潮汐強(qiáng)度及海床地形等因素的影響，不同區(qū)域鹽水濃度不同[2]，同時(shí)，不同季節(jié)由于徑流存在差別，鹽水上溯距離也不同[3]。由于影響因素眾多，因此鹽淡水混合是極其復(fù)雜的科學(xué)問題。

目前，國內(nèi)外關(guān)于河口鹽淡水混合的研究成果主要聚焦在機(jī)理和鹽水入侵及其工程應(yīng)用上。在鹽淡水混合機(jī)理方面，徑潮流的此消彼長(zhǎng)作用是影響鹽淡水平面分布、垂向混合的主要因素[4-5]，此外地形、工程及其氣象條件等因素對(duì)鹽淡水混合也有較大影響[6]；鹽水入侵研究表明，除外界因素外，鹽淡水混合后導(dǎo)致的水體密度的變化也會(huì)對(duì)水流流態(tài)產(chǎn)生影響[7]，不同河口的影響程度不同；鹽水入侵工程應(yīng)用目前關(guān)注較多的是河口水源地飲用水安全問題，例如針對(duì)長(zhǎng)江口青草沙水庫，眾多學(xué)者在分析研究海平面上升、南水北調(diào)、極端天氣、人類活動(dòng)、河勢(shì)變化等對(duì)青草沙水庫取水產(chǎn)生影響的基礎(chǔ)上，提出了利用三峽水庫補(bǔ)淡壓咸、北支建閘等工程措施[8-9]，并建立了長(zhǎng)江口鹽水入侵預(yù)報(bào)系統(tǒng)[10]，為水源地取水安全提供預(yù)警、預(yù)報(bào)。

本文在前人研究成果的基礎(chǔ)上，分析椒江河口鹽水入侵特征。椒江位于錢塘江南側(cè)約150 km 處（位置見圖1），是我國典型的強(qiáng)潮河口，河口口門處海門站多年平均潮差為4.0 m，最大潮差可達(dá)6.3 m[11]。通過本文研究，一方面可以明晰強(qiáng)潮河口的鹽水入侵特征，另一方面可以為相關(guān)科學(xué)研究以及工程建設(shè)提供技術(shù)支撐和依據(jù)。

圖1 模型范圍與網(wǎng)格示意圖Fig.1 Schematic diagram of model domain and grid

1 鹽水入侵?jǐn)?shù)學(xué)模型構(gòu)建

1.1 模型范圍

本文重點(diǎn)研究椒江河口徑流與潮流相互作用下的鹽水上溯規(guī)律，因此首先建立平面二維鹽水入侵?jǐn)?shù)學(xué)模型。數(shù)值模擬軟件采用南京水利科學(xué)研究院自主研發(fā)的CJK3D-WEM 數(shù)值模擬系統(tǒng)[12]，該軟件適用于江河湖泊、河口海岸等涉水工程中的水動(dòng)力、泥沙、水質(zhì)、溫排、溢油的模擬預(yù)測(cè)。笛卡爾坐標(biāo)系下的水流以及鹽水運(yùn)動(dòng)二維方程表達(dá)如下：

式中：H為水深，單位：m；z為水位，單位：m；u、v為速度分量，單位：m/s；t為時(shí)間，單位：s；f為科氏力，單位：s-1；g為重力加速度，單位：m/s2；Nx、Ny為x、y向水流紊動(dòng)粘性系數(shù)，單位：m2/s；ρ0為水體密度，單位：kg/m3；τsx、τsy為x、y向的風(fēng)切應(yīng)力，τbx、τby為x、y向的底切應(yīng)力，單位：kg/（m·s2）；?為鹽度，單位：‰；Kx、Ky為鹽水?dāng)U散系數(shù)，單位：m2/s。

采用三角形網(wǎng)格有限體積法求解式（1），對(duì)流項(xiàng)數(shù)值通量采用Roe 格式的近似Riemann 解，紊動(dòng)項(xiàng)采用單元交界面的平均值來計(jì)算通過該界面紊動(dòng)粘性項(xiàng)的數(shù)值通量，底坡項(xiàng)采用斜底方法處理，方程的具體求解過程見文獻(xiàn)[13]。

數(shù)學(xué)模型構(gòu)建以椒靈江干流和臺(tái)州灣為核心水域，北側(cè)邊界至寧波象山縣，南側(cè)邊界至樂清灣附近，東側(cè)邊界至椒江口門以東85 km 處的開闊水域，上游邊界包含椒靈江一級(jí)支流永安溪和始豐溪（見圖1）。模型上游邊界采用流量控制，外海邊界采用潮位控制。模型網(wǎng)格總數(shù)為50 572 個(gè)，椒靈江水域網(wǎng)格加密，最小網(wǎng)格邊長(zhǎng)約為23.5 m，最大網(wǎng)格邊長(zhǎng)為3 508 m，椒靈江干流河床地形采用2016 年實(shí)測(cè)地形，臺(tái)州灣及外海水深數(shù)據(jù)采用最新海圖水深數(shù)據(jù)，統(tǒng)一到85 高程，坐標(biāo)系統(tǒng)采用2000 國家大地坐標(biāo)系。數(shù)學(xué)模型計(jì)算時(shí)間步長(zhǎng)取1 s；科氏力f=2ωsinφ，其中ω為地球自轉(zhuǎn)角速度，φ為地理緯度，取28.7°；糙率n采用附加糙率方式進(jìn)行處理，計(jì)算公式為：

式中：H為水深。水流紊動(dòng)粘性系數(shù)取0.01HU*（U*為摩阻流速）；動(dòng)邊界水深取0.02 m；鹽水?dāng)U散系數(shù)取30HU*。

1.2 模型驗(yàn)證

數(shù)學(xué)模型驗(yàn)證包括水動(dòng)力和鹽度兩部分。采用2016 年椒靈江河道同步水文測(cè)驗(yàn)資料，包括3 個(gè)潮位站和5條垂線（見圖2），驗(yàn)證時(shí)間為2016年12月8日00 時(shí)（北京時(shí)，下同）—16 日00 時(shí)，包括大潮、中潮、小潮的連續(xù)過程。圖3—4為潮位以及潮流流速流向驗(yàn)證。從圖中可以看出，數(shù)學(xué)模型的潮位計(jì)算值與實(shí)測(cè)值吻合較好，西門、西岙、海門站高低潮位偏差基本在0.1 m 以內(nèi)，大潮、中潮、小潮的潮流流速及流向與實(shí)測(cè)值基本吻合，流速偏差基本在10%以內(nèi)，流向偏差小于15°，數(shù)學(xué)模型水動(dòng)力驗(yàn)證總體滿足相關(guān)規(guī)程要求。圖5 為各垂線鹽度驗(yàn)證，由于實(shí)測(cè)資料鹽度采樣時(shí)間跨度（3 h）較大，因此鹽度的數(shù)學(xué)模型計(jì)算值無法與實(shí)際鹽度逐時(shí)對(duì)比，但總體來看，兩者大小及過程均吻合良好，鹽度自下游向上游逐漸減小，5#點(diǎn)最高鹽度值在20‰左右，3#點(diǎn)最高約為12‰，1#點(diǎn)為2.5‰左右。從水動(dòng)力與鹽度的驗(yàn)證結(jié)果來看，本文建立的數(shù)學(xué)模型能夠較好地模擬椒江鹽水運(yùn)動(dòng)特征。

圖2 潮位站及潮流垂線位置示意圖Fig.2 Schematic diagram of tide level station and tide vertical line position

圖3 水位過程驗(yàn)證Fig.3 Tide level verification

圖4 潮流過程驗(yàn)證Fig.4 Tidal current verification

圖5 鹽水入侵過程驗(yàn)證Fig.5 Salinity intrusion verification

2 椒江河口鹽水入侵特征研究

2.1 計(jì)算條件

椒江為典型的山溪性河流，流量年內(nèi)分配不均，月度差異明顯，根據(jù)保證率為50%的流域來水條件分析（見圖6），11月平均流量最小，約為23 m3/s，6 月最大，約為684 m3/s，全年平均流量約為160 m3/s，洪枯季流量差異可達(dá)30 倍以上。為反映不同徑流下的鹽水運(yùn)動(dòng)特征，上游徑流分別取10 m3/s、100 m3/s、200 m3/s、500 m3/s、1 000 m3/s和3 000 m3/s，外海潮位采用1.2 節(jié)驗(yàn)證大潮、中潮、小潮潮位過程的結(jié)果，海門站大潮平均潮差約為5.2 m，小潮平均潮差約為2.8 m，根據(jù)2020 年海門站潮差累積頻率統(tǒng)計(jì)進(jìn)行計(jì)算，得到大潮累積頻率約為5%，小潮約為87%（見圖6）。

圖6 2001年椒江徑流月度分配（a）與2020年海門站潮差累積頻率統(tǒng)計(jì)（b）Fig.6 Monthly distribution of Jiaojiang River runoff in 2001（a）and statistics of tidal range cumulative frequency of Haimen station in 2020（b）

2.2 鹽水分布特征

鹽淡水混合是徑潮流共同作用的結(jié)果。圖7與表1 給出了不同徑流、不同潮差條件下鹽水的最大上溯距離分布（指距離頭門島的距離）。從中可以看到，大潮時(shí)鹽水上溯距離明顯大于小潮，隨著徑流的增大，鹽水逐漸向下游運(yùn)動(dòng)，當(dāng)徑流為10 m3/s時(shí)，0.1‰鹽水可以上溯至始豐溪與永安溪的交匯處三江村附近，5‰鹽水可以上溯至八仙巖附近；當(dāng)徑流為200 m3/s 時(shí)，0.1‰鹽水僅能上溯至西岙附近，距離口外頭門島（見圖2中C1處）約60 km，5‰鹽水上溯至海門上游13 km 處；當(dāng)徑流為500 m3/s 時(shí)，0.1‰鹽水僅能上溯至石仙婦上游約6 km 處，距離口外頭門島約40 km，5‰鹽水上溯至海門上游4 km處；當(dāng)徑流為3 000 m3/s時(shí)，海門以上河段的鹽水濃度均小于0.1‰。

表1 不同條件下的鹽水最大上溯距離與鹽水鋒面比降Tab.1 Maximum upstream tracking distance of Brine and the gradient of salt water front under different conditions

圖7 不同徑流條件下的鹽水最大上溯距離Fig.7 Maximum upstream tracking distance of brine under different runoff conditions

在椒江干流沿程布置間距為5～5.5 km 的17個(gè)采樣點(diǎn)（位置見圖2）用于分析不同潮差、不同徑流條件下的椒江鹽水鋒面分布，結(jié)果見圖8。由圖可知，無論是小潮還是大潮，低流量條件下由于鹽水上溯距離遠(yuǎn)，鹽水鋒面坡度較緩，當(dāng)徑流量為10 m3/s時(shí)，小潮的鹽水濃度比降約為0.52‰/km，大潮約為0.51‰/km；隨著徑流量的增大，鹽水濃度比降逐漸增大（見表1），當(dāng)徑流量為500 m3/s，小潮的鹽水濃度比降約為0.93‰/km，大潮約為0.73‰/km；當(dāng)徑流量為1 000 m3/s，小潮的鹽水濃度比降約為0.94‰/km，大潮約為0.76‰/km。綜上可見，隨著徑流量的增大，鹽水達(dá)到最大上溯距離后，鹽水鋒面比降變化梯度呈減小趨勢(shì)。

圖8 不同徑流潮差條件下的鹽水鋒面分布Fig.8 Salt water front distribution under different runoff and tidal ranges

2.3 徑流對(duì)鹽水上溯距離變化梯度影響分析

在特定的徑流條件下，鹽水在漲落潮的作用下最大上溯位置不斷發(fā)生變化，漲憩時(shí)鹽水上溯距離最遠(yuǎn)，落憩時(shí)最近。圖9 統(tǒng)計(jì)了椒江河口鹽水在大潮和小潮不同潮型作用下的鹽水上溯距離變化梯度與徑流變化的關(guān)系，即：

圖9 鹽水上溯距離變化與徑流的關(guān)系Fig.9 Relationship between saltwater upstream tracking distance change and runoff

式中：ε 為鹽水上溯距離隨徑流的變化梯度，單位：m/（m3·s）；?Q為不同徑流Q1、Q2的差值，單位：m3/s；?d為徑流Q1、Q2下鹽水上溯距離d1、d2之間的差值，單位：m。

由圖9可以看到，無論在漲憩還是落憩時(shí)刻，隨著徑流量的增大，鹽水上溯距離呈減速降低趨勢(shì)，漲憩時(shí)刻鹽水上溯距離變化梯度略大于落憩，小潮鹽水上溯距離變化梯度大于大潮。在10～100 m3/s的流量區(qū)間，漲憩時(shí)刻大潮鹽水上溯距離變化梯度為222 m/（m3·s），小潮為238 m/（m3·s），落憩時(shí)刻兩者數(shù)值分別為207 m/（m3·s）和233 m/（m3·s），而在200～500 m3/s流量區(qū)間，鹽水上溯距離變化梯度急劇降低至15～35 m/（m3·s），這說明在徑流抑制鹽水上溯的過程中，鹽水鋒面距離河口越遠(yuǎn)，其受徑流的影響越大，反之，受徑流影響相對(duì)較弱。

2.4 潮差對(duì)鹽水上溯變化區(qū)間影響分析

根據(jù)數(shù)學(xué)模型計(jì)算結(jié)果，統(tǒng)計(jì)了在特定徑流條件下的小潮、中潮和大潮（對(duì)應(yīng)潮型25 h 內(nèi)的最大潮差分別為3.2 m、5.1 m、6.2 m）鹽水上溯距離變化區(qū)間（漲憩鹽水上溯距離-落憩鹽水上溯距離），結(jié)果見圖10。由圖可知，隨著潮差的增大，鹽水上溯變化區(qū)間總體呈線性增長(zhǎng)趨勢(shì)，隨著徑流量的增大，鹽水上溯變化區(qū)間逐漸減小。當(dāng)徑流為10 m3/s時(shí)，小潮鹽水上溯變化區(qū)間約為12 km，大潮約為21 km，而徑流為500 m3/s 時(shí)，大小潮的鹽水上溯變化區(qū)間明顯下降，小潮約為8 km，大潮約為14 km，說明鹽水鋒面距離河口越遠(yuǎn)，潮差對(duì)鹽水上溯變化區(qū)間影響越大，反之，鹽水鋒面距離河口越近，潮差對(duì)鹽水上溯變化區(qū)間影響越小。

圖10 潮差對(duì)鹽水上溯區(qū)間影響Fig.10 Influence of tidal range on salt water upstream tracking region

2.5 椒江河口鹽水上溯距離函數(shù)關(guān)系式

自然條件下，鹽水上溯距離主要與徑流和潮差有關(guān)。根據(jù)表1中徑流與鹽水上溯距離的關(guān)系分析可知，在潮差不變的條件下，鹽水上溯距離與徑流的自然對(duì)數(shù)（lnQ）具有一定線性關(guān)系，因此用lnQ代表徑流因子，同時(shí)在徑流不變的條件下，鹽水上溯距離與潮差的線性關(guān)系較好（見圖10）。由以上分析可以構(gòu)建椒江河口鹽水最大上溯距離與徑流、潮差的函數(shù)關(guān)系式：

式中：d為鹽水上溯距離（即距離頭門島的距離，單位：km）；Q為徑流，單位：m3/s；?h為潮差，單位：m；a、b、c為系數(shù)，由于鹽水上溯距離與徑流呈反比，與潮差呈正比，因此a＜0，b＞0。

根據(jù)數(shù)學(xué)模型計(jì)算不同徑流、不同潮差下漲憩和落憩各18 組鹽水上溯的最大距離，利用式（4）進(jìn)行擬合，求解a、b、c。經(jīng)擬合計(jì)算，漲憩時(shí)a= -9.95，b= 2.18，c= 89.67；落憩時(shí)a= -8.08，b= 0.34，c= 75.96。圖11為鹽水上溯距離函數(shù)關(guān)系式預(yù)測(cè)值與數(shù)學(xué)模擬值的對(duì)比結(jié)果，經(jīng)統(tǒng)計(jì)，漲憩時(shí)刻兩者偏差為2%～3%，落憩時(shí)刻為5%～6%。漲憩和落憩時(shí)刻數(shù)學(xué)模擬值與預(yù)測(cè)值的相關(guān)性高達(dá)0.99，說明本文建立的椒江河口鹽水上溯距離函數(shù)關(guān)系式預(yù)測(cè)精度相對(duì)較高，可為椒江鹽水相關(guān)研究提供借鑒。

圖11 鹽水上溯距離函數(shù)關(guān)系式預(yù)測(cè)值與數(shù)學(xué)模擬值對(duì)比Fig.11 Comparison between predicted value and digital analog value of salt water upstream tracking distance function

3 結(jié)語

以強(qiáng)潮河口——椒江河口為研究對(duì)象，通過數(shù)學(xué)模型研究了不同徑流、潮差條件下的鹽水入侵特征。主要結(jié)論如下：

①椒江河口是典型的強(qiáng)潮河口，當(dāng)徑流為10 m3/s 時(shí)，0.1‰鹽水可以上溯至始豐溪與永安溪的交匯處三江村附近，5‰鹽水可以上溯至八仙巖附近，當(dāng)徑流為500 m3/s 時(shí)，0.1‰鹽水僅能上溯至石仙婦上游約6 km 處，5‰鹽水上溯至海門上游4 km 處，當(dāng)徑流為3 000 m3/s時(shí)，海門以上椒江河段的鹽水濃度均小于0.1‰。

②在徑流變化條件下，鹽水上溯距離變化梯度漲憩大于落憩，小潮大于大潮，隨著徑流量的增大，鹽水上溯距離變化梯度逐漸增大；在潮差變化條件下，鹽水上溯距離變化區(qū)間呈線性增長(zhǎng)趨勢(shì)，隨著徑流量的增大，同等潮差下的鹽水上溯距離變化區(qū)間逐漸減小。

③通過數(shù)學(xué)模擬結(jié)果擬合出了椒江河口鹽水最大上溯距離與徑流、潮差的函數(shù)關(guān)系式，漲憩時(shí)刻關(guān)系式預(yù)測(cè)值偏差為2%～3%，落憩時(shí)刻為5%～6%，預(yù)測(cè)值精度較高，可為椒江鹽水相關(guān)研究提供借鑒。