汾河蘭村—柴村段河流水位抬升對潛水流場影響的數值模擬研究

王瑞芬 季晨龍 李海明

(1.山西省水利水電科學研究院，山西 太原 030002；2.天津科技大學濱海地下水利用與保護研究室，天津 300457；3.天津市海洋資源與化學重點實驗室，天津 300457)

由于經濟社會的快速發展和人口的急劇增長，汾河流域生態環境受到嚴重破壞[1]。山西省先后組織了4次大規模治理，基本實現了汾河干流全年不斷流，流域地下水位持續回升，但流域自身的產流功能仍未恢復[2]。汾河流域歷史上擁有豐富的濕地資源，在涵養水源、降解污染、防洪蓄水等方面起到重要作用[3]，但在汾河早期的治理中，原本屬于河道的大量灘涂、濕地被開發利用，濕地資源不斷減少[4]。為此，山西省委省政府實施了汾河流域生態修復中游核心區干流蓄水工程，該工程主要通過河道蓄水向堤外濕地補水，蓄水期間使流域內濕地水位抬升，蓄水水位保持在灘面1～2m。

流域內濕地水位的抬升使得在濕地補給地下水的過程中面臨各種技術與科學問題，例如濕地水位的升降容易引起土地鹽堿化、沼澤化等問題。目前針對蓄水工程對下游地下水的影響多從滲流量[5]、環境流量[6]角度展開，也有研究者運用不同的方法，對汾河流域的地下水位進行了動態預測[7-10]，但對濕地水位抬升影響下的地下水流場特征研究較為罕見。

數值模擬也已成為研究河水與地下水交換的重要手段[11-14]。其中，對河流與含水層系統模型的正確概化以及對模型參數的準確估算是影響研究結果的關鍵因素[15-16]。本文運用數值模擬的方法，研究了汾河蘭村—柴村段在河流水位瞬時抬升0.5m、1.0m、2.0m后的潛水流場特征，可為山西省汾河流域濕地水資源管理和修復提供相關參考。

1 研究區水文地質條件

研究區位于太原市西張盆地，具體位置見圖1。西、北兩面被群山環抱，海拔多在780m以上，面積約為66km2。蘭村—柴村段的河流屬平原性河流，長度約10km，寬約50～300m，隨著河流流程的增加，含水層富水性增大。地下水主要接受汾河干流段地表水的滲漏補給以及大氣降水入滲補給，地下水位低于河水位，河水通過入滲向河兩岸含水層側向滲透補給地下水。潛水的地下水徑流主要受地形、含水層巖性、補給來源等控制。總體上由北部向南部，由東西兩側向中部流動運移。

圖1 研究區位置

研究區地下水為松散巖類孔隙水，含水層多為全新統砂卵礫石層，在0～160m埋深范圍內分為潛水含水層和弱透水層，包氣帶厚度約10～22m，巖性多為粉土、亞砂土。潛水含水層頂板高程為620～824m，厚度約30～150m，巖性多為砂卵石。弱透水層厚度約10～80m，底板埋深40～160m，巖性多為黏土。

潛水的排泄方式主要有蒸發、越流排泄及含水層的側向徑流3種形式。

2 潛流帶水動力特征數值模擬

2.1 水文地質概念模型的建立

2.1.1 邊界條件

根據地形地貌、地層結構及區域流場特征，潛水含水層南邊界與東邊界中下部概化為定流量邊界，東邊界上部、西邊界以及北邊界概化為隔水邊界，弱透水層均概化為隔水邊界，見圖2。

圖2 潛水含水層邊界條件

2.1.2 含水層空間概化

本文只關注與河流水力聯系密切的潛水含水層，根據已有的地質圖、剖面圖、地質報告等資料將含水層概化為非均質各向同性潛水含水層，地下水流以水平運動為主。研究區含水層三維結構與剖面概化見圖3。

圖3 研究區含水層三維結構與剖面概化

2.2 地下水流數值模型的建立

2.2.1 空間剖分

對模擬研究區進行網格剖分，共剖分為20000個網格(100行×100列×2層)，其中活動單元格為8368個，垂向上劃分為兩層，見圖4。

圖4 網格剖分示意圖

2.2.2 模型識別應力期選取

根據現有資料，模型以2012年1—12月為校正期，2013年1—12月為驗證期，以月為應力期，每個應力期為1個時間步長。

2.2.3 源匯項處理

本文模型中潛水含水層地下水補給項主要包括降水入滲補給、地表水滲漏補給，排泄項主要包括蒸發、徑流排泄。

降水入滲補給量根據降水量乘以降水入滲補給系數得到；地表水滲漏補給根據水文站徑流量測量資料，按一定滲漏比計算得到；潛水蒸發根據已有資料，潛水含水層現狀地下水位埋深在臨界埋深(5m)以下的區域不考慮潛水蒸發，在臨界埋深以上的沖積、洪積交界洼地區潛水蒸發系數值采用0.040，沖積平原區采用0.013；邊界徑流排泄量采用達西定律計算得到。

2.3 模型的校正與驗證

根據研究區校正期和驗證期潛水流場模擬值與實測值的比較，兩者的擬合程度較高，說明所建立的水文地質概念模型和數值模型基本合理，該模型能較好地模擬研究區地下水動態變化過程。

2.4 情景假設

本文以2014年1月地下水位為初始條件，預測了河流水位瞬時抬升0.5m、1.0m、2.0m時引起潛水含水層地下水位的變化情況，并且在每個河流抬升幅度下選取5日、10日、30日、60日、90日來舉例說明。

預測期模型主要考慮的補給項有降水入滲補給、地表水滲漏補給，排泄項主要有蒸發、徑流排泄。其中降雨量采用2004—2014年降雨量的平均值作為降水入滲補給輸入模型，其他參數與2013年保持不變。

3 結果與討論

研究區初始潛水流向為北部由北向南；中部受河流影響，流向變為由西北向東南；南部為由東西兩側向中間流動。北部等水位線密集，潛水流速快，水力坡度約為0.0065，流速約為0.220m/d。中部與南部等水位線稀疏，流速與北部相比較緩慢，水力坡度約為0.0009，流速約為0.041m/d。研究區內河流水位高于潛水水位，河流與潛水的補給關系為河流補給潛水。河流水位瞬時抬升0.5m、1.0m、2.0m后，在5日、10日時對潛水流場整體影響不大，但河流處的等水位線出現微小變化，河水補給兩側潛水，河流對潛水的補給作用增強，見圖5、圖6。

圖5 5日潛水模擬流場

圖6 10日潛水模擬流場

隨著時間的推移，河流對潛水的側滲補給不斷加強，河流水位瞬時抬升0.5m、1.0m、2.0m時，河流附近的潛水流向逐漸變為從河流向東西兩岸流動，特別是當河流水位抬升1.0m、2.0m時變化最為明顯，見圖7、圖8。

圖7 30日潛水模擬流場

圖8 60日潛水模擬流場

河流水位瞬時抬升0.5m、1.0m、2.0m之后，經過90日，研究區潛水流向在由北向南的趨勢上，由于河流的側滲作用，河流附近的潛水流向明顯變為由河流向東西兩岸補給，見圖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日時影響距離已經超出研究區邊界。河流水位瞬時抬升后，左岸90日時部分區域的影響距離也已超出研究區邊界。

圖12 河流水位抬升2.0m影響距離

根據圖10～圖12，得出河流抬升不同幅度時，經過不同時間，其相應的影響距離，見表1。多數情況下，河流水位抬升幅度越高，相同時刻下的潛水影響距離越大。

表1 不同時間不同河水位抬升幅度潛水水位影響距離 單位：m

4 結 論

a.研究區初始潛水流向為北部由北向南；中部受河流影響，流向變為由西北向東南；南部為由東西兩側向中間流動。北部等水位線密集，潛水流速快，水力坡度約為0.0065，流速約為0.220m/d。中部與南部等水位線稀疏，流速與北部相比較緩慢，水力坡度約為0.0009，流速約為0.041m/d。河流水位高于潛水水位，河流與潛水的補給關系為河流補給潛水。

b.隨著時間的推移以及河流水位抬升幅度的增加，河流對潛水的側滲補給不斷加強，使河流附近的潛水流向最先開始改變，經過90日之后，潛水流向在由北向南的趨勢上，河流附近的潛水流向已明顯變為由河流向東西兩岸補給。

c.多數情況，河流水位抬升相同幅度時，隨著時間的推移，影響距離逐漸增加；河流水位抬升幅度越高，相同時刻下的潛水影響距離越大。河流水位瞬時抬升0.5m、1.0m、2.0m時，河流左岸60日影響距離分別為1742～2345m、1541～2613m、1809～3015m；90日時部分區域已超出研究區范圍。