王瑞芬 季晨龍 李海明
(1.山西省水利水電科學研究院,山西 太原 030002;2.天津科技大學濱海地下水利用與保護研究室,天津 300457;3.天津市海洋資源與化學重點實驗室,天津 300457)
由于經(jīng)濟社會的快速發(fā)展和人口的急劇增長,汾河流域生態(tài)環(huán)境受到嚴重破壞[1]。山西省先后組織了4次大規(guī)模治理,基本實現(xiàn)了汾河干流全年不斷流,流域地下水位持續(xù)回升,但流域自身的產(chǎn)流功能仍未恢復[2]。汾河流域歷史上擁有豐富的濕地資源,在涵養(yǎng)水源、降解污染、防洪蓄水等方面起到重要作用[3],但在汾河早期的治理中,原本屬于河道的大量灘涂、濕地被開發(fā)利用,濕地資源不斷減少[4]。為此,山西省委省政府實施了汾河流域生態(tài)修復中游核心區(qū)干流蓄水工程,該工程主要通過河道蓄水向堤外濕地補水,蓄水期間使流域內(nèi)濕地水位抬升,蓄水水位保持在灘面1~2m。
流域內(nèi)濕地水位的抬升使得在濕地補給地下水的過程中面臨各種技術(shù)與科學問題,例如濕地水位的升降容易引起土地鹽堿化、沼澤化等問題。目前針對蓄水工程對下游地下水的影響多從滲流量[5]、環(huán)境流量[6]角度展開,也有研究者運用不同的方法,對汾河流域的地下水位進行了動態(tài)預測[7-10],但對濕地水位抬升影響下的地下水流場特征研究較為罕見。
數(shù)值模擬也已成為研究河水與地下水交換的重要手段[11-14]。其中,對河流與含水層系統(tǒng)模型的正確概化以及對模型參數(shù)的準確估算是影響研究結(jié)果的關(guān)鍵因素[15-16]。本文運用數(shù)值模擬的方法,研究了汾河蘭村—柴村段在河流水位瞬時抬升0.5m、1.0m、2.0m后的潛水流場特征,可為山西省汾河流域濕地水資源管理和修復提供相關(guān)參考。
研究區(qū)位于太原市西張盆地,具體位置見圖1。西、北兩面被群山環(huán)抱,海拔多在780m以上,面積約為66km2。蘭村—柴村段的河流屬平原性河流,長度約10km,寬約50~300m,隨著河流流程的增加,含水層富水性增大。地下水主要接受汾河干流段地表水的滲漏補給以及大氣降水入滲補給,地下水位低于河水位,河水通過入滲向河兩岸含水層側(cè)向滲透補給地下水。潛水的地下水徑流主要受地形、含水層巖性、補給來源等控制??傮w上由北部向南部,由東西兩側(cè)向中部流動運移。
圖1 研究區(qū)位置
研究區(qū)地下水為松散巖類孔隙水,含水層多為全新統(tǒng)砂卵礫石層,在0~160m埋深范圍內(nèi)分為潛水含水層和弱透水層,包氣帶厚度約10~22m,巖性多為粉土、亞砂土。潛水含水層頂板高程為620~824m,厚度約30~150m,巖性多為砂卵石。弱透水層厚度約10~80m,底板埋深40~160m,巖性多為黏土。
潛水的排泄方式主要有蒸發(fā)、越流排泄及含水層的側(cè)向徑流3種形式。
2.1.1 邊界條件
根據(jù)地形地貌、地層結(jié)構(gòu)及區(qū)域流場特征,潛水含水層南邊界與東邊界中下部概化為定流量邊界,東邊界上部、西邊界以及北邊界概化為隔水邊界,弱透水層均概化為隔水邊界,見圖2。
圖2 潛水含水層邊界條件
2.1.2 含水層空間概化
本文只關(guān)注與河流水力聯(lián)系密切的潛水含水層,根據(jù)已有的地質(zhì)圖、剖面圖、地質(zhì)報告等資料將含水層概化為非均質(zhì)各向同性潛水含水層,地下水流以水平運動為主。研究區(qū)含水層三維結(jié)構(gòu)與剖面概化見圖3。
圖3 研究區(qū)含水層三維結(jié)構(gòu)與剖面概化
2.2.1 空間剖分
對模擬研究區(qū)進行網(wǎng)格剖分,共剖分為20000個網(wǎng)格(100行×100列×2層),其中活動單元格為8368個,垂向上劃分為兩層,見圖4。
圖4 網(wǎng)格剖分示意圖
2.2.2 模型識別應力期選取
根據(jù)現(xiàn)有資料,模型以2012年1—12月為校正期,2013年1—12月為驗證期,以月為應力期,每個應力期為1個時間步長。
2.2.3 源匯項處理
本文模型中潛水含水層地下水補給項主要包括降水入滲補給、地表水滲漏補給,排泄項主要包括蒸發(fā)、徑流排泄。
降水入滲補給量根據(jù)降水量乘以降水入滲補給系數(shù)得到;地表水滲漏補給根據(jù)水文站徑流量測量資料,按一定滲漏比計算得到;潛水蒸發(fā)根據(jù)已有資料,潛水含水層現(xiàn)狀地下水位埋深在臨界埋深(5m)以下的區(qū)域不考慮潛水蒸發(fā),在臨界埋深以上的沖積、洪積交界洼地區(qū)潛水蒸發(fā)系數(shù)值采用0.040,沖積平原區(qū)采用0.013;邊界徑流排泄量采用達西定律計算得到。
根據(jù)研究區(qū)校正期和驗證期潛水流場模擬值與實測值的比較,兩者的擬合程度較高,說明所建立的水文地質(zhì)概念模型和數(shù)值模型基本合理,該模型能較好地模擬研究區(qū)地下水動態(tài)變化過程。
本文以2014年1月地下水位為初始條件,預測了河流水位瞬時抬升0.5m、1.0m、2.0m時引起潛水含水層地下水位的變化情況,并且在每個河流抬升幅度下選取5日、10日、30日、60日、90日來舉例說明。
預測期模型主要考慮的補給項有降水入滲補給、地表水滲漏補給,排泄項主要有蒸發(fā)、徑流排泄。其中降雨量采用2004—2014年降雨量的平均值作為降水入滲補給輸入模型,其他參數(shù)與2013年保持不變。
研究區(qū)初始潛水流向為北部由北向南;中部受河流影響,流向變?yōu)橛晌鞅毕驏|南;南部為由東西兩側(cè)向中間流動。北部等水位線密集,潛水流速快,水力坡度約為0.0065,流速約為0.220m/d。中部與南部等水位線稀疏,流速與北部相比較緩慢,水力坡度約為0.0009,流速約為0.041m/d。研究區(qū)內(nèi)河流水位高于潛水水位,河流與潛水的補給關(guān)系為河流補給潛水。河流水位瞬時抬升0.5m、1.0m、2.0m后,在5日、10日時對潛水流場整體影響不大,但河流處的等水位線出現(xiàn)微小變化,河水補給兩側(cè)潛水,河流對潛水的補給作用增強,見圖5、圖6。
圖5 5日潛水模擬流場
圖6 10日潛水模擬流場
隨著時間的推移,河流對潛水的側(cè)滲補給不斷加強,河流水位瞬時抬升0.5m、1.0m、2.0m時,河流附近的潛水流向逐漸變?yōu)閺暮恿飨驏|西兩岸流動,特別是當河流水位抬升1.0m、2.0m時變化最為明顯,見圖7、圖8。
圖7 30日潛水模擬流場
圖8 60日潛水模擬流場
河流水位瞬時抬升0.5m、1.0m、2.0m之后,經(jīng)過90日,研究區(qū)潛水流向在由北向南的趨勢上,由于河流的側(cè)滲作用,河流附近的潛水流向明顯變?yōu)橛珊恿飨驏|西兩岸補給,見圖9。北部水力坡度與潛水流速變幅較小,中部、南部水力坡度和流速與初始流場相比均有降低,水力坡度為0.0004~0.0006,流速為0.017~0.027m/d。
圖9 90日潛水模擬流場
圖10 河流水位抬升0.5m影響距離
圖11 河流水位抬升1.0m影響距離
將相同時刻下河流水位抬升影響后的潛水流場與初始流場相比較,繪制出河流水位抬升不同幅度時對潛水的影響距離,見圖10~圖12。由圖10~圖12可以看出,河流水位抬升相同幅度下,隨著時間的推移,影響距離逐漸增加,并且由于河流右岸空間有限,在河流水位抬升0.5m 60日、河流水位抬升1.0m以及2.0m 30日時影響距離已經(jīng)超出研究區(qū)邊界。河流水位瞬時抬升后,左岸90日時部分區(qū)域的影響距離也已超出研究區(qū)邊界。
圖12 河流水位抬升2.0m影響距離
根據(jù)圖10~圖12,得出河流抬升不同幅度時,經(jīng)過不同時間,其相應的影響距離,見表1。多數(shù)情況下,河流水位抬升幅度越高,相同時刻下的潛水影響距離越大。
表1 不同時間不同河水位抬升幅度潛水水位影響距離 單位:m
a.研究區(qū)初始潛水流向為北部由北向南;中部受河流影響,流向變?yōu)橛晌鞅毕驏|南;南部為由東西兩側(cè)向中間流動。北部等水位線密集,潛水流速快,水力坡度約為0.0065,流速約為0.220m/d。中部與南部等水位線稀疏,流速與北部相比較緩慢,水力坡度約為0.0009,流速約為0.041m/d。河流水位高于潛水水位,河流與潛水的補給關(guān)系為河流補給潛水。
b.隨著時間的推移以及河流水位抬升幅度的增加,河流對潛水的側(cè)滲補給不斷加強,使河流附近的潛水流向最先開始改變,經(jīng)過90日之后,潛水流向在由北向南的趨勢上,河流附近的潛水流向已明顯變?yōu)橛珊恿飨驏|西兩岸補給。
c.多數(shù)情況,河流水位抬升相同幅度時,隨著時間的推移,影響距離逐漸增加;河流水位抬升幅度越高,相同時刻下的潛水影響距離越大。河流水位瞬時抬升0.5m、1.0m、2.0m時,河流左岸60日影響距離分別為1742~2345m、1541~2613m、1809~3015m;90日時部分區(qū)域已超出研究區(qū)范圍。