變化環(huán)境下萬泉河河口區(qū)鹽分擴散的模擬分析

海來以波，施浩然，劉曉慶，劉 蛟，方明儀，張藍月

(1.西華大學能源與動力工程學院，四川 成都 610039；2.西南科技大學環(huán)境與資源學院，四川 綿陽 621010)

河口區(qū)域在內(nèi)河和外海因素交互影響下形成特殊的水生態(tài)環(huán)境，是洄游魚類產(chǎn)卵、育幼和索餌的重要場所[1-2]。隨著上游水資源開發(fā)帶來的內(nèi)河來水條件變化[3-4]，以及海平面上升等帶來的河口潮位波動[5]，導致外海潮波帶來的鹽分在河口區(qū)的擴散發(fā)生變化，水體鹽度的改變對河口區(qū)水生態(tài)環(huán)境產(chǎn)生了重要的影響。張洪亮等[6]對萊州灣的研究說明鹽度升高引起了浮游植物生物量降低、浮游生物的種類減少，同時對海洋生物產(chǎn)卵場造成了不利影響。除了徑流和潮位2個主導因素外，河口區(qū)鹽度尤其是橫斷面上的分布同樣受到河床地形的影響[7]。林一楠等[8]通過對錢塘江河口沿程7個測站的實測數(shù)據(jù)分析，揭示了鹽水主要沿深槽上溯、凹岸鹽度及其變幅大于凸岸的客觀規(guī)律。

河口區(qū)的鹽分擴散是由內(nèi)河來流、外海潮位和河床地形三者共同作用下形成的水動力條件決定的。因此，基于水動力學模型的研究成為河口區(qū)咸潮入侵防治以及生態(tài)環(huán)境修復的主要手段。Duong等[9]采用Mike21模型，通過一維和二維的耦合模型模擬了越南湄公河三角洲一條主要支流在未來不同情景下(如海平面上升、上游來流量變化等)的海水入侵變化。Uncles等[10]對強潮作用下英格蘭南部Lynher河口徑流和海水入侵相互作用的季節(jié)變化進行了研究，結(jié)果表明海水入侵主要受潮汐和徑流的控制。趙高磊等[11]以咸潮上溯和沖淡區(qū)為研究對象，分析河口上游修建水庫后，在豐水期、平水期、枯水期徑流改變的情況下，北門江河口地區(qū)咸潮上溯和沖淡區(qū)變化規(guī)律。Xu等[12]基于Mike21分析了河道排放條件、潮汐范圍和海水入侵之間的關(guān)系，構(gòu)建了長江口的咸水入侵函數(shù)。陳浩等[13]基于Mike21 FM 建立了研究區(qū)域海水入侵模型，對辛集擋潮閘不同下泄流量時原始河道和清淤后河道的海水入侵進行模擬研究，以論證所估算生態(tài)基流的合理性。

萬泉河為海南省第三大河，近年來由于發(fā)展需要，河道水庫的修建和運行改變了河口區(qū)原有的水文條件[4,14]，同時中國近海區(qū)域海平面的上升也影響了萬泉河河口潮位變化。徑流和潮位使得萬泉河河口區(qū)鹽分擴散的變化，直接影響到了對鹽度極為敏感的國家二級保護魚類花鰻鱺的生境條件[15-16]，同時也對河口區(qū)域內(nèi)的紅樹林生態(tài)林以及農(nóng)田生態(tài)系統(tǒng)造成影響，沿海農(nóng)田土壤鹽漬化問題更加凸顯[17- 18]。為此，本文基于Mike21模型，開展徑流和潮位變化情景下，不同水文時期河道內(nèi)鹽度分布模擬，分析萬泉河特殊的河網(wǎng)結(jié)構(gòu)下，不同季節(jié)中徑流和潮位變化對海水潮波上行距離、影響范圍以及強度等方面的具體作用，闡明不同條件下萬泉河河口區(qū)的鹽分擴散變化規(guī)律，為河口區(qū)的水環(huán)境保護、水生態(tài)建設(shè)提供依據(jù)。

1 研究區(qū)域概況

萬泉河為海南省第三大河流，發(fā)源于瓊中黎族苗族自治縣五指山，經(jīng)萬寧市至瓊海市，在博鰲港匯入南海，干流全長157 km，平均坡降1.12‰，總落差586 m。萬泉河流域位于109°37′E～110°38′E和18°46′N～19°31′N，整體形狀呈帶狀，平均寬度15 km，面積3 693 km2。本文對外海潮波上行中鹽分在內(nèi)河河道擴散進行研究，因此根據(jù)多年來萬泉河河口區(qū)潮波影響范圍的潮區(qū)界確定本次研究的范圍為入海口至嘉積鎮(zhèn)大坡村，河流長度約為16 km。研究區(qū)河道主要流經(jīng)嘉積鎮(zhèn)和博鰲鎮(zhèn)，屬于國家重點開發(fā)區(qū)，是海南省生態(tài)空間格局中以旅游業(yè)為主題的東部核心區(qū)，也是亞洲論壇永久性會址所在地(圖1)。

圖1 研究區(qū)域

萬泉河河口區(qū)的河網(wǎng)結(jié)構(gòu)復雜，尤其是在河流入?？诟浇膷u嶼眾多，干流被主要的沙坡島和東嶼島分為北側(cè)、中側(cè)和南側(cè)三支，與河口聯(lián)成一體，形成較大水域，而該水域與外海由橫亙的海沙與河沙沖積而成的玉帶灘隔斷，僅在北端有一個小的入?？?，為此在萬泉河河口形成特殊的葫蘆形港灣。

萬泉河流域河口的潮汐類型屬混合潮，每一太陰日內(nèi)有時發(fā)生2次高潮和2次低潮，有時只出現(xiàn)1次高潮和1次低潮，潮汐日不等現(xiàn)象顯著。根據(jù)河口處附近港北潮位站的監(jiān)測數(shù)據(jù)，按國家85高程基面，河口的多年平均高潮位和低潮位分別為 0.96、0.21 m，潮汐差0.75 m。位于研究區(qū)上游約5 km的嘉積水文站為萬泉河下游的控制性水文站，監(jiān)測的多年平均流量為147 m3/s。萬泉河河口區(qū)上游的天然來流量年內(nèi)分配極不均勻，汛期6—11月占了全年徑流量的80%左右，9—11月的主汛期占60%左右，但是隨著上游水資源的開發(fā)，經(jīng)上游水庫調(diào)節(jié)過后，河口區(qū)的年平均來流量和年內(nèi)分布差異均有所減少。

2 數(shù)值模擬

2.1 Mike21模型

由于本次研究河口區(qū)域水深較淺，加上縱向水流對鹽度分布的擾動，鹽度分層不明顯，因此本次研究根據(jù)收集的數(shù)據(jù)，采用丹麥DHI公司開發(fā)的Mike21模型對萬泉河河口區(qū)在2017年3個典型水文月中日尺度上的鹽度分布開展模擬。Mike21主要包含二維水動力模塊、波浪模塊、簡單對流擴散模塊、水生態(tài)模塊以及泥沙運移模塊等，在河口、海灣和海岸地區(qū)的水動力及其相關(guān)現(xiàn)象的平面二維仿真模擬得到廣泛使用[11,13,19-20]。Mike21水動力模型遵循納維耶-斯托克斯(Navier-Stokes)方程，并服從Boussinesq假設(shè)，即流體低速流動中，密度的變化不考慮由壓強變化所引起的，僅考慮溫度對密度的影響，同時服從靜水壓力。Mike21的控制方程組如下：

(1)

(2)

Mike21模型采用有限體積法，對研究區(qū)進行三角網(wǎng)格化離散處理。本次模型劃分為7 690個網(wǎng)格。網(wǎng)格最小角度29°，最大角度59.9°，共計4 388個計算節(jié)點，計算的最小時間步長取1.5 s，最大時間步長取200 s，CFL(Courant-Friedrichs-Lewy)數(shù)為0.8。

2.2 模型數(shù)據(jù)

按流量大小，對萬泉河河口區(qū)2017年5、2、4月3個不同水文月的鹽度分布開展模擬，模型上邊界條件的驅(qū)動數(shù)據(jù)為模擬期內(nèi)嘉積水文站的逐日流量數(shù)據(jù)；下邊界條件的驅(qū)動數(shù)據(jù)為模擬期內(nèi)港北潮位站的逐日潮位數(shù)據(jù)。降水、蒸發(fā)、風速風向等氣象數(shù)據(jù)來自于國家氣象科學數(shù)據(jù)中心(http://data.cma.cn/)的瓊海氣象站實測數(shù)據(jù)。

模擬河道的河床組成主要為中粗砂，床面不平整，采砂形成的局部深坑散亂分布，灘地較寬闊生長有雜草和低矮灌木，兩側(cè)邊岸一般為茂密灌木和椰子、檳榔等高大喬木林帶，河道兩岸為較為光滑的堤防工程，根據(jù)河道現(xiàn)狀，結(jié)合河床地形分布，研究區(qū)域的糙率取為0.025～0.030。

研究河段內(nèi)無實測的水文數(shù)據(jù)，且考慮不同徑流和潮位條件下的鹽度分布，而鹽度分布是鹽分在水動力學作用下的結(jié)果，因此直接根據(jù)汀州水質(zhì)監(jiān)測站的監(jiān)測斷面在各月份的實測鹽度數(shù)據(jù)開展模型驗證。

2.3 模擬情景

首先根據(jù)2017年的實測數(shù)據(jù)，對代表性水文月的鹽度模擬結(jié)果進行率定和驗證，并將其作為參照模型。在此基礎(chǔ)上，考慮上游來流變化和河口潮位變化2種情景，通過改變參照組的輸入條件，設(shè)置對比模型。具體情景設(shè)置如下。

a)參照組模型。萬泉河河口區(qū)多年平均流量為147 m3/s，該流量對河道的塑造以及穩(wěn)定生態(tài)環(huán)境的形成具有重要的意義。同時根據(jù)水文的季節(jié)性變化，結(jié)合上游的水庫調(diào)節(jié)作用，選擇了2017年5、2、4月，開展模擬與驗證。3個月的平均流量分別為178、148、106 m3/s，平均高潮位分別為0.96、0.93、0.96 m，平均低潮位分別為0.15、0.12、0.18 m。

b)流量對比組模型?？紤]研究區(qū)域上游在水資源開發(fā)中的取水增加后，在參照組的基礎(chǔ)上，分別假設(shè)河口區(qū)的來流減少10%和20%，而潮位不變的2種情景，對河口區(qū)鹽分擴散進行模擬，分析內(nèi)河來流對萬泉河河口區(qū)鹽分擴散的影響。

c)潮位對比組模型。基于中國近海近50年海平面變化的周期性、趨勢性等規(guī)律，預測研究河口區(qū)的海平面可能具有較大幅度的上升[21]。因此，在參照組的基礎(chǔ)上，假設(shè)潮位增加0.1、0.3 m，上游來流量不變的2種情景，對河口區(qū)鹽分擴散進行模擬，量化分析外海潮位對萬泉河河口區(qū)鹽分擴散的影響。

3 模型驗證與結(jié)果分析

3.1 模型驗證

根據(jù)3個月份模擬的平均鹽度與實測數(shù)據(jù)的對比(表1)看出，3個月的鹽度誤差不大，最大的5月為-13.9%，4月的誤差僅為3.3%。該結(jié)果說明構(gòu)建Mike21模型能夠較好地反映出萬泉河河口區(qū)在不同徑流和潮位作用下的鹽分擴散情況。

表1 汀州斷面鹽度的模擬與實測結(jié)果對比

3.2 鹽分沿河擴散距離變化

本文參照相關(guān)研究以河道內(nèi)鹽度濃度為0.5‰位置為外海鹽分上溯擴散頂點對上溯距離進行統(tǒng)計[11]，得到了不同情景下萬泉河河口區(qū)咸潮沿河上溯的影響范圍，結(jié)果見表2。總體上來看，上游來流量的減少以及外海潮位的上升均會對鹽分沿河擴散距離產(chǎn)生正向影響，且來流變化的影響表現(xiàn)更加明顯。在不同月份，當流量發(fā)生同等幅度的減小時，流量偏豐月份的咸潮上溯距離變化更大。當潮位增加0.1 m時，就會對來流量為106 m3/s的4月份鹽分上溯距離產(chǎn)生影響，而在來流量為148 m3/s的2月時，潮位需要增加0.3 m才會對鹽分上溯距離產(chǎn)生比較明顯的影響；當來流量增大至178 m3/s的5月時，假設(shè)情景下潮位的增加對鹽分上溯距離基本沒有影響。通過不同情景下，鹽分上溯距離的變化可以看出，當上游來流量偏枯的時候，外海潮汐變化的作用更加明顯；而當上游來流量偏豐的時候，內(nèi)河來流量變化的作用更加明顯。

表2 萬泉河河口區(qū)不同情景下鹽分沿河擴散距離統(tǒng)計

3.3 河道內(nèi)鹽度分布變化

2017年4月萬泉河河口區(qū)以及4種變化情景下河道內(nèi)的鹽度分布見圖2，可以看出不同情境下，河道的鹽度分布變化規(guī)律基本一致，在3條支流上都表現(xiàn)為鹽度在0.8 km范圍內(nèi)快速地從4‰降低到0.5‰以下。這主要是受4月份上游來流量較少的影響，內(nèi)河徑流對入?？诟浇}分擴散的影響不大。因此，在徑流量減小的2種情景下，1‰以上的鹽度分布區(qū)域基本沒有變化，只是上溯距離分別增加了0.26、0.28 km。而潮波通過入?？谶M入內(nèi)河后，首先較大的水域起到了緩沖作用，降低了潮汐作用；同時萬泉河河口區(qū)在入海口附近河床平緩，而3條支流的河床地形往上游抬升明顯，這進一步降低了潮汐對鹽分的上行擴散作用，因此在潮位升高的2種情景下，1‰以上的鹽度分布區(qū)域也基本沒有變化，只是上溯距離分別增加了0.21、0.38 km。

圖2 2017年4月萬泉河河口區(qū)不同情景鹽度分布

相對于4月來水量，5月萬泉河的流量偏豐，其不同情景下河道內(nèi)的鹽度分布見圖3，可以看出不同來流減少以及潮位升高后鹽分上溯頂點(0.5‰)的上移變化，原有流量178 m3/s減少10%和20%后，外海鹽分的影響范圍分別增加了0.25、0.52 km；而隨著潮位增加0.1、0.3 m，咸水上溯距離無明顯變化。從鹽度的空間分布變化上來看，總體上不同情景的分布規(guī)律較為類似，表現(xiàn)為在狹小的入?？诟浇鼌^(qū)域鹽度最高，均高于4‰，在繼續(xù)沿河擴散中隨著內(nèi)河水域的擴大，鹽度比較均勻地往上游擴散后逐步減小，然而繼續(xù)向上游擴散進入分叉河道后，受彎道水動力學塑造的河床地形影響，在三支分叉河道中表現(xiàn)出不同的規(guī)律。首先在北側(cè)支流河岸凸點下游，河道的流量大，流速大，該區(qū)域的水深較其他地區(qū)淺，鹽度擴散快，因此鹽度較低，為0.1‰～0.5‰，而在上游，受彎曲河道頂托作用，沿河流速減少，鹽分在此處累積，呈現(xiàn)出比下游更高的鹽度分布。同樣在南側(cè)支流，東嶼島與玉帶灘對中間河道的束窄，鹽分在潮汐作用下通過束窄區(qū)進入上游后，受到增強的內(nèi)河徑流影響，鹽分在此處累積，鹽度達到4‰左右。在經(jīng)過分叉河道，進入汀州斷面上游的主河道之后，河道內(nèi)的鹽度擴散又呈現(xiàn)出比較均勻的擴散，直至上溯頂點。

圖3 2017年5月萬泉河河口區(qū)不同情景鹽度分布

4 討論

通過對5種情景下萬泉河河口鹽度分布的對比模擬分析，可以看出，外海鹽分在河口區(qū)的擴散影響范圍受到了潮汐增加和來流減小的正向影響，這與目前基本的研究成果是一致的[10-11,22]。但是在萬泉河口，不同程度的流量減小都會使得鹽分擴散范圍增加，但是潮位升高的影響還受到了內(nèi)河徑流的制約。當萬泉河上游來流為106 m3/s，潮位增加0.1 m就能使得鹽分擴散距離上移，而來流量增大到148 m3/s，潮位需要增加0.3 m才能使得鹽分擴散距離上移。這說明徑流和潮汐對鹽分擴散范圍的影響是相互影響、相互制約的。

研究發(fā)現(xiàn)，萬泉河河口區(qū)特殊的水網(wǎng)形狀以及河床地形等對徑流和潮汐相互作用下形成的水動力條件具有決定性的作用，這在很大程度上決定了河口區(qū)的鹽度分布。這與Parsa等[23]得出的河口的鹽度入侵在很大程度上取決于河口的形態(tài)、地形、水動力和潮汐特性的結(jié)論相一致。萬泉河在出?？谛纬傻乃?qū)Τ毕鸬搅司彌_作用，使得大部分鹽分聚集在此區(qū)域，鹽度分布迅速降低；加之島嶼的阻隔與河床抬升的影響，這都削弱了潮汐對鹽分的擴散影響作用。這就可能造成了模擬結(jié)果中潮位變化情景下的鹽度分布改變?nèi)跤趶搅髯兓?。同時由于河道的彎曲和束窄變化，在局部形成特殊的水流特征，出現(xiàn)鹽分在局部累積的特殊分布情況。

本文從對河道環(huán)境等具有重要意義的流量出發(fā)，考慮一定的水文變化，在不同的來流和潮位變化情景下，對萬泉河河口在2017年3個典型水文月份上的鹽度分布情況進行了模擬分析，闡明了在特殊的水系結(jié)構(gòu)及河道地形條件下，徑流和潮位變化對萬泉河河口鹽度分布的影響作用，可為河口區(qū)的水生態(tài)環(huán)境建設(shè)提供科學依據(jù)。