淺析采用三維數值分析方法計算水庫地下水壅高

裴 穎

（長安大學環境科學與工程學院 西安 710021）

1 概述

水庫蓄水后，庫區周圍地下水相應壅高而接近或高于地面，會產生一些不良現象，如農田鹽漬化、沼澤化及建筑物條件惡化、礦坑涌水等。水庫蓄水浸沒問題是水庫工程勘察與評價的重要內容之一，是確定庫區征地和移民遷建水位的重要依據。嚴重的水庫浸沒問題影響到水庫正常蓄水位的選擇，甚至影響到壩址的選擇。浸沒評價結論是水庫設計、建設非常重要的因素和先決條件，如果水庫蓄水后的地下水埋深值小于浸沒地下水埋深臨界值，就會產生浸沒。

江巷水庫位于江淮分水嶺北部、滁州市定遠縣連江鎮境內，其主要任務是供水和灌溉，兼顧防洪等綜合利用。水庫總庫容為1.30億m3，水庫正常蓄水位43.0m。工程規模為大（2）型。洪水標準為100年一遇設計，2000年一遇校核，設計、校核洪水位分別為44.42m和45.50m。

江巷水庫為典型的平原地區水庫，庫周地勢較為低緩，易產生浸沒問題，可研設計階段需開展浸沒問題研究，通過一系列勘探、試驗、觀測等方法進行浸沒評價，判定水庫蓄水后是否會產生浸沒問題。作為重要判定指標的庫岸地下水位壅高值可通過地下水動力學方法計算，主要方法有解析法、三維數值分析法等。本文主要通過三維數值分析法來計算水庫地下水位壅高值，得出庫岸地下水埋深，為水庫的浸沒評價提供依據。

2 三維滲流模型計算地下水壅高方法

2.1 數學模型

假定庫區為非均質各向同性，三維地下水流非穩定運動的數學模型可表示為：

式中：Ω為滲流區域，H為含水層的水位標高，K為滲透系數，μ為給水度；ε為含水層的源匯項，H0為含水層的初始水位分布，Γ1為滲流區域的二類邊界，包括承壓含水層底部隔水邊界和滲流區域的側向流量或隔水邊界，n為邊界面的法向方向，Kn為邊界面法向方向的滲透系數，q（x，y，z，t)定義為二類邊界的單位面積流量，流入為正、流出為負，隔水邊界為0。

2.2 有限元離散

對式（1）進行有限元離散，得到整個滲流區域的代數方程：

式中：[G]為總體滲透矩陣（傳導矩陣），[P]為貯水矩陣，｛H｝為未知節點水頭列陣，｛F｝為已知的右端項。

采用“隱式差分格式”，有：

整理后，得：

2.3 水文地質參數確定

庫區地層較為簡單，表層均為第四系上更新統（Q3a）l地層覆蓋，以重粉質壤土、粉質粘土為主，下伏基巖為白堊系紅層，地表無出露。在庫區郭小圩、孟嚴、東方曹、齊崗、小廟陳布置了5個地質剖面（圖1），

圖1 庫區地質剖面布置圖

庫區地下水類型主要為基巖裂隙水和孔隙水兩大類。含水層位主要為下伏厚層紅層碎屑巖層頂部風化帶裂隙水和上覆粘性土孔隙潛水。地下水賦存條件很差，水量較為貧乏。

（1）滲透系數

根據地質資料表明，庫區均為第四系土層覆蓋，層厚15～30m，以重粉質壤土、粉質粘土為主，該土層滲透系數一般都在n×10-6～n×10-7cm/s之間。下伏基巖為泥質、鈣質粉砂巖和砂質泥巖，亦不易形成滲透通道，不具備形成地下大規模強滲流的地質條件，可以當做隔水邊界處理。另外，通過對庫區的地質調查發現，庫區內所有的沖溝都是由第四系土層組成，為微透水或極微透水土層，構成庫盆的地層都是第四系上更新統（Q3al）的粘性土，地下不存在強滲透含水層。滲透系數取值參考注水試驗結果取值（表1）。

表1 滲透系數及壓力傳導系數取值表

（2）給水度

庫區的巖性主要為重粉質壤土、粉質粘土，給水度取值為5%。

（3）降水入滲補給系數

降水入滲補給系數α是指降水滲入量與降水總量的比值，α值的大小取決于地表土層的巖性和土層結構、地形坡度、植被覆蓋以及降水量的大小和降水形式等，它是一個無量綱系數，其值變化于0～1之間，庫區的年均降雨量約為934mm，主要巖性為重粉質壤土、粉質粘土，降水入滲補給系數取值為0.1。

（4）潛水蒸發系數

潛水蒸發系數主要年水面蒸發量、巖性和地下水位埋深有關。庫區的多年年水面蒸發量為798mm，地下水位埋深為0.5～1.5m，主要巖性為重粉質壤土、粉質粘土，蒸發系數取值為0.18。潛水蒸發量主要與潛水位的埋深、包氣帶巖性、地表植被和氣候等因素相關。根據前人的研究成果，一般認為水位埋深大于5m的地區潛水蒸發很小。庫區地下水位埋深一般小于5m，因此要考慮蒸發影響。

2.4 水文地質概念模型

水文地質概念模型是在綜合分析地下水系統的基礎上，對評價區地質、含水層實際的邊界條件、內部結構、滲透性質、水力特征和補給排泄等水文地質條件進行科學的綜合、歸納和加工，從而對一個復雜的水文地質實體進行概化，便于進行數學或者物理模擬。因此，建立水文地質概念模型主要應該考慮如下幾個方面：概化后的模型應該具備反映庫區水文地質原型的功能，概化后的各類邊界條件應符合庫區地下水流場特征，概化后的模型邊界應該盡量利用自然邊界，人為邊界性質的確定應從不利因素考慮等。本次模型采用46m等高線作為模型邊界，按照水文地質勘察鉆孔所監測的地下水位插值獲得邊界水位（表2和圖2）。

表2 庫區觀測孔地下水位統計表

圖2 通過觀測孔的地下水位插值得到的庫區地下水位等值線圖

模型四周為水頭邊界，在43m以下高程的也為水頭邊界，頂部為降雨入滲和蒸發邊界，模型底部基巖為隔水邊界，建立相應水文地質概念模型。

2.5 區域離散

計算區域以項目所在地中心位置為坐標原點，正北方向為y軸正向，正東方向為x軸正向，垂直向上為z軸正向，垂向上考慮1層，將庫區域離散為51001個節點，57860個單元。

2.6 邊界和初始條件

邊界條件：庫區為一個相對獨立的水文地質單元，地下水由周邊向庫區流動，外邊界為水頭邊界，通過實測地下水位插值獲得；內邊界分別為高程低于常水位42.5m和正常蓄水位43m的區域，上部邊界為降雨入滲和蒸發邊界，底部為隔水邊界。初始條件：考慮穩定流計算。

3 庫岸剖面地下水壅高計算

在庫區選取了5個地質剖面計算水庫蓄水后的地下水位變化，本文列出其中的兩個地質剖面的計算結果（圖3，其他剖面情況類似）。從圖中可以看出，不同剖面地下水位存在一定差別，但總體上地下水位浸潤線比較平緩。

圖3 庫岸剖面地下水位等值線圖

4 計算結果分析

根據蓄水位為43m時的地下水位的分布情況，庫區周邊地下水位略高于庫水位，地下水對庫水位進行補給，其水位差不大，水力坡度較小。但大壩下游地下水等值線密集，水位降幅較大，水力梯度較大，表明大壩防滲性能較好。本次同時也計算了水庫水位為常水位42.5m時的地下水位的分布情況，其地下水位分布規律與水庫蓄水位為43m時基本一致。

三維數值分析法計算結果表明，與初始地下水位相比，蓄水后地下水壅高為0.09～1.53m。當水庫蓄水位為42.5m時，庫區周邊區域地下水位高于常水位0.32～0.56m，平均值為0.41m。當水庫蓄水位為43m時，庫區周邊區域地下水位高于常水位0.48～0.66m，平均值為0.52m。這是在庫水位長期穩定在42.5m和43m時的地下水位變化，實際庫水位一直在變化，因此，穩定庫水位條件下計算的地下水位是最終的理想結果，比實際地下水位偏高。

從計算結果可以看出庫岸地下水位的變化規律和趨勢：水庫運行多年后，從水庫向四周，地下水位逐漸升高。當庫水位上升時，與庫水邊距離越近，地下水位上升越快；距離越遠，地下水位上升越慢。地下水位上升（壅高）較高的位置一般在庫水邊，越遠離庫水邊上升越少。雖然不同區域土層參數有所差異，影響范圍和水位上升高度也有差異，但地下水位上升趨勢總體一致，地下水位浸潤線比較平緩，庫岸地下水位隨庫水位上升需要很長時間才能達到穩定。

5 結語

庫岸地下水位壅高值計算主要采用地下水動力學方法，較常用的有解析法和三維數值分析法。兩種方法各有其適用條件和特點。與解析法相比，三維數值分析法更適用于滲流場較復雜的浸沒區，通過對可能產生浸沒的區域內地下水位情況進行模擬，可以有效解決預測的浸沒范圍與實際浸沒范圍相比偏大的問題。對于規模等別較高的水庫，可分別采用解析法和三維數值分析法進行浸沒計算，對比、驗證、分析計算結果，對更準確、合理地提出浸沒評價結論是大有裨益的■