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

海來以波，施浩然，劉曉慶，劉 蛟，方明儀，張藍(lán)月

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

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

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

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

1 研究區(qū)域概況

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

圖1 研究區(qū)域

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

萬泉河流域河口的潮汐類型屬混合潮，每一太陰日內(nèi)有時(shí)發(fā)生2次高潮和2次低潮，有時(shí)只出現(xiàn)1次高潮和1次低潮，潮汐日不等現(xiàn)象顯著。根據(jù)河口處附近港北潮位站的監(jiān)測(cè)數(shù)據(jù)，按國家85高程基面，河口的多年平均高潮位和低潮位分別為 0.96、0.21 m，潮汐差0.75 m。位于研究區(qū)上游約5 km的嘉積水文站為萬泉河下游的控制性水文站，監(jiān)測(cè)的多年平均流量為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ū)域水深較淺，加上縱向水流對(duì)鹽度分布的擾動(dòng)，鹽度分層不明顯，因此本次研究根據(jù)收集的數(shù)據(jù)，采用丹麥DHI公司開發(fā)的Mike21模型對(duì)萬泉河河口區(qū)在2017年3個(gè)典型水文月中日尺度上的鹽度分布開展模擬。Mike21主要包含二維水動(dòng)力模塊、波浪模塊、簡(jiǎn)單對(duì)流擴(kuò)散模塊、水生態(tài)模塊以及泥沙運(yùn)移模塊等，在河口、海灣和海岸地區(qū)的水動(dòng)力及其相關(guān)現(xiàn)象的平面二維仿真模擬得到廣泛使用[11,13,19-20]。Mike21水動(dòng)力模型遵循納維耶-斯托克斯(Navier-Stokes)方程，并服從Boussinesq假設(shè)，即流體低速流動(dòng)中，密度的變化不考慮由壓強(qiáng)變化所引起的，僅考慮溫度對(duì)密度的影響，同時(shí)服從靜水壓力。Mike21的控制方程組如下：

(1)

(2)

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

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

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

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

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

2.3 模擬情景

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

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

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

c)潮位對(duì)比組模型?；谥袊＝?0年海平面變化的周期性、趨勢(shì)性等規(guī)律，預(yù)測(cè)研究河口區(qū)的海平面可能具有較大幅度的上升[21]。因此，在參照組的基礎(chǔ)上，假設(shè)潮位增加0.1、0.3 m，上游來流量不變的2種情景，對(duì)河口區(qū)鹽分?jǐn)U散進(jìn)行模擬，量化分析外海潮位對(duì)萬泉河河口區(qū)鹽分?jǐn)U散的影響。

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

3.1 模型驗(yàn)證

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

表1 汀州斷面鹽度的模擬與實(shí)測(cè)結(jié)果對(duì)比

3.2 鹽分沿河擴(kuò)散距離變化

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

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

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)河徑流對(duì)入?？诟浇}分?jǐn)U散的影響不大。因此，在徑流量減小的2種情景下，1‰以上的鹽度分布區(qū)域基本沒有變化，只是上溯距離分別增加了0.26、0.28 km。而潮波通過入海口進(jìn)入內(nèi)河后，首先較大的水域起到了緩沖作用，降低了潮汐作用；同時(shí)萬泉河河口區(qū)在入海口附近河床平緩，而3條支流的河床地形往上游抬升明顯，這進(jìn)一步降低了潮汐對(duì)鹽分的上行擴(kuò)散作用，因此在潮位升高的2種情景下，1‰以上的鹽度分布區(qū)域也基本沒有變化，只是上溯距離分別增加了0.21、0.38 km。

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

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

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

4 討論

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

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

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