平原水庫滲流場模擬及滲漏量計算
——以斗門水庫試驗段為例

李文溢， 楊阿敏， 周維博

(1.長安大學 環境科學與工程學院， 陜西 西安 710054；2.長安大學 旱區地下水文與生態效應教育部重點實驗室， 陜西 西安 710054)

1 研究背景

平原水庫是在平原區天然湖泊、洼地、河道等地形的基礎上，通過圈筑圍堤而建成的水庫[1]。在我國西北干旱及半干旱地區，由于受到人類活動和經濟發展等多方面的影響，區域性的水環境問題、生態環境問題日益突出[2]，嚴重影響了地區社會發展及“一帶一路”建設，而平原水庫在解決這些問題和促進區域經濟發展等方面起著重要作用。但平原水庫的運行也往往會對區域地下水、地表水、土壤環境等帶來一定的負面影響，產生庫水富營養化、泥沙淤積、浸沒、土壤鹽漬化等一系列問題[3]，這對水庫的運行效率和城市供水安全造成了嚴重的隱患。

針對平原水庫運行引起的環境問題，國內外學者展開了大量研究。李榮榮等[4]采用遙感影像法分析了水庫下游土壤鹽分指數變化，結果表明水庫蓄水使得下游土壤鹽漬化程度不斷加重；Winter[5]在充分考慮介質各向異性、地下水流系統、滲透系數等因素的情況下，運用垂向二維穩定流模型研究了水庫和地下水之間的相互作用；韓菲等[6]系統分析了國內外湖泊及水庫富營養化模型的優缺點，結果表明富營養化模型應當著重考慮水動力、化學、生物3方面因素的影響。

目前關于平原水庫的研究多集中在庫水水質評價及周邊土壤生態環境等方面，對水庫滲漏及周邊地下水流場演變等方面的研究較少，而這正是威脅平原水庫安全運行的常見問題[7-9]。在天然條件下，地下水系統會保持一定的均衡狀態，而當水庫蓄水之后，庫水的滲漏補給作用改變了區域地下水系統的收支特征，打破了原有的地下水均衡動態[10]；另外，庫底沉積物的分布、滲透性特征、地下水位變化會對水庫的滲漏速率產生影響[11]。因此在充分考慮水庫沉積物滲透性、地下水位特征的基礎上，分析計算水庫滲漏量，探究水庫蓄水后區域地下水流場的變化就顯得尤為重要。

本文以陜西省西咸新區斗門水庫試驗段為例，采用數值法和解析法相結合，分析庫區滲漏量及水庫蓄水對區域地下水流場的影響，以期為庫底防滲、水庫調蓄、地下水環境保護及斗門水庫的進一步建設提供相關科學依據。

2 研究區概況

斗門水庫在陜西省原漢代昆明池遺址的基礎上修建而成，位于西咸新區灃東新城境內(原西安市長安區)，北接魚斗路，西靠新韋斗路，距離灃河約2.6 km。斗門水庫是引漢濟渭工程的重要一環，在輸配水工程中起著重要作用，是一座集調蓄、防洪、改善生態環境等多項功能于一體的水資源綜合利用平原水庫[12]。斗門水庫試驗段是斗門水庫工程建設的一部分，于2017年3月正式蓄水運行，庫區面積0.53 km2，其中水面面積約0.47 km2，庫底高程396.8 m，正常蓄水位400.47 m，蓄水深度3.67 m，蓄水量155×104m3，地理位置如圖1所示。

圖1 研究區地理位置圖

2.1 地形地貌

斗門水庫試驗段所處地貌單元為秦嶺北麓近山前沖洪積平原，庫區內渭河一級階地大面積分布，階面高程為396～410 m，整體上東南高、西北低，地形開闊，高程起伏較為平緩。

2.2 地層巖性

庫區地表50 m以內地層主要形成于第四系，可分為以下幾種類型：(1)上更新統(Qp3)風積層，其巖性主要為風積黃土(Q32eol)，層厚約5～12 m，外表呈淡灰黃色，垂直節理及大孔隙發育，質地較為疏松，分布于庫區西南部。(2)全新統(Qh)沖積層，分為上下兩部。下部沖積層(Q41al)分上下兩層：上層厚約10～40 m，主要為壤土、粉土，夾少量中砂，于庫區南部、西部(渭河一、二級階地)廣泛分布；下層巖性主要為粉質黏土，局部夾薄砂層，夾層厚度0.3～0.5 m。而上部沖積層(Q42al)巖性主要為壤土，層厚約5～10 m，受到人類活動影響，頂部人工堆積零星分布(約0.5 m，主要為水泥塊、磚塊等建筑垃圾)，于庫區東北部(渭河一級階地)分布。

2.3 水文地質條件

庫區地下水類型為孔隙潛水，埋深在7.5～14.2 m之間，主要賦存于第四系黏性土層(夾部分少量砂層)中，大致由東南流向西北，含水層埋藏于山前洪積扇、階地及漫灘。受到斗門水庫建設影響，周邊農田面積及人為灌溉有所減少，因此降水入滲補給和人工開采分別為斗門水庫區域地下水的主要補給和排泄途徑。

3 庫區滲漏量計算

3.1 水庫滲漏模型

在水庫的建設過程中，庫壩體內側面做了較好的防滲處理，因此可以忽略內側面的滲漏，僅考慮庫底區域，根據庫底沉積物巖性將水庫劃分為不同的滲漏分區進行分析計算。庫水對地下水入滲補給的滲漏模型如圖2所示，可概化為：(1)庫水在水頭差的作用下通過庫底沉積物垂向入滲到含水層中；(2)在地下水流場的作用下，發生側向滲流。根據達西定律[13-14]，采用公式(1)所示的離散化模型計算庫區滲漏量：

(1)

式中：Q為庫區總滲漏量，即庫水補給地下水的總量,m3/d；Qi為分區滲漏量，即庫水在各個分區對地下水的補給量,m3/d；Ai為各分區庫底的水平面積,m2；ΔHi為庫水位與所對應的地下水監測井水位之差,m；Li為監測井與所對應的分區形心的距離，m；Ki為各分區沉積物的垂向滲透系數，m/d；n為分區個數。

圖2 水庫滲漏模型示意圖

3.2 水庫滲漏分區

根據斗門水庫試驗段工程地質勘察資料，按照庫底沉積物巖性，將庫區劃分為3個區域，滲漏分區及監測井位示意圖如圖3所示。其中Ⅰ-Ⅲ為不同的滲漏分區，GW1-GW6為沿庫北岸自西向東布設的6口地下水位監測井。

圖3 庫區滲漏分區圖

考慮到水庫較深，不適宜進行野外原位試驗測定沉積物滲透系數。根據土工試驗標準[15]，在庫區內共取樣12組，采用室內變水頭滲透試驗法[16]測定了庫底沉積物的滲透系數。變水頭滲透試驗法原理如圖4所示：將野外所取的原狀土樣作為滲流通道，將變水頭管的供水作為水源來模擬地下水的一維滲流過程；通過多次調節(5～6次)不同的水頭高度進行重復試驗；記錄多次試驗的時間和水頭來計算分析試驗結果。

圖4 變水頭滲透試驗法示意圖

根據試驗數據，采用公式(2)計算庫底沉積物的垂向滲透系數：

(2)

圖4及公式(2)中：a為變水頭管垂直于水流方向的橫斷面面積,cm2；L為變水頭管兩端過水斷面之間的距離(滲徑)，在本試驗裝置中等于土樣高度,cm；KV為沉積物的垂向滲透系數,m/d；t1、t2分別為同一次試驗中初始、最后讀取變水頭管水頭的時間,s；H1、H2分別為前后兩次讀取的變水頭管的水頭,m。

通過野外取樣及室內試驗求得了各個分區庫底沉積物的垂向滲透系數。各滲漏分區面積、庫底沉積物巖性及垂向滲透系數如表1所示。

由表1可以看出：庫底沉積物滲透系數量級為10-3m/d，其中庫區東北部(Ⅱ區)沉積物滲透系數相對較大。

表1 滲漏分區相關參數

3.3 水庫滲漏計算

結合表1中的相關參數，根據庫區2018年3-7月的地下水位監測資料(由于庫區施工限制，4月份監測資料缺失)，求得了庫區滲漏量及日滲漏強度[17]，結果如表2所示。

表2 庫區滲漏量及滲漏強度

由表2可以看出：庫區日平均滲漏量為19.37 m3/d，由此可得年滲漏量約為7 070.05 m3/a，約占水庫總庫容(155×104m3)的4.6‰，庫區平均日滲漏強度約為0.041 mm/d，年滲漏強度約為14.97 mm/a。從時間上看：3-7月庫區滲漏量變化呈下降趨勢，分析原因是由于3-7月隨著降雨量的增大，庫區地下水位有所抬升，而由于水庫的運行調節，庫水位基本無變化，因此使得區域水力坡度減小，庫水對地下水的滲漏補給減小；從空間上看：庫區東北部(II區)的滲漏強度較大，分析其原因是II區庫底沉積物巖性以淤泥質土為主，且夾有部分中砂，因此造成了較大的滲漏損失，存在較大的滲漏風險。在水庫運行管理中，應當著重加強庫區東北側的地下水位監測，做好防滲措施。

4 地下水流數值模擬

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

4.1.1 模擬范圍 由于斗門水庫試驗段面積較小，大尺度的區域水流模擬不能精確反映斗門水庫試驗段蓄水對周邊地下水流場的影響。因此采用等流量交換法[18]，概化模擬區范圍。如圖5所示：模擬區北至常七隊，南到楊家莊村，距離約1.42 km，西接西安龍門中學，東至袁其寨一帶，距離約1.85 km，總面積約1.83 km2。模擬計算的目標含水層為第四紀孔隙潛水含水層，厚度約55～60 m，含水層上邊界為潛水面，下邊界為隔水黏土層。

4.1.2 水力特征及邊界條件概化 模擬區內含水層地下水主要表現為水平運動，僅在河道、庫岸周邊等區域存在垂向運動，整體地下水流速較小，符合層流運動的特征，適用于Darcy定律；滲透系數、給水度等水文地質參數在含水層不同點存在一定差異，為非均質性；地下水系統的補給、徑流、排泄等要素隨時間和空間有所變化，為非穩定流。綜合分析模擬區的水文地質條件及地下水流的補、徑、排條件，將模型邊界條件概述如下：將黏性土隔水底板概化為潛水底部隔水邊界；潛水面接受大氣降水補給、庫水入滲補給、蒸發排泄等垂向入滲補給，為不穩定水面；將東北及西南側邊界概化為零流量邊界，東南側及西北側為二類流量邊界。

4.2 地下水流系統數學模型

根據斗門水庫試驗段模擬區的水文地質概念模型，可建立如式(3)所示的數學模型來描述模擬區含水層地下水的二維水流運動：

(3)

式中：x，y為坐標變量；D為模擬區范圍；h為水位,m；h0為初始水位,m；K為滲透系數,m/d；W為垂向補給強度,m/d；μ為給水度；Zb為含水層底板標高，m;q為二類邊界上單位面積流入(流出)水量，m3/(d·m2);Γ1、Γ2分別為一、二類邊界；n為二類邊界外法線方向。

4.3 模型設計

4.3.1 模型網格剖分 根據實際水文地質條件和計算需求，采用不等距剖分的方法對所建立的模型進行了網格剖分，將模型剖分成了若干行和列，在庫岸邊進行了一定程度的網格加密。如圖6所示：共剖分132行、160列，21 120個單元格，其中活動單元格15 789個，非活動單元格5 331個。

4.3.2 模型頂板和底板高程 通過DEM數據提取獲取模擬區頂板高程數據，通過收集模擬區鉆孔資料及工程地質勘察資料獲取模擬區底板高程數據，采用克里金法進行插值，生成模擬區頂、底板高程如圖7所示。

4.3.3 模型源匯項的處理 根據相關研究[19-20]，潛水埋深超過5 m，其蒸發可忽略不計，而模擬區潛水埋深在9～14 m之間，因此在模型源匯項的處理中不考慮潛水蒸發，主要考慮以下源匯項：(1)降雨入滲補給及蒸發。根據模擬區氣象站實測降雨蒸發資料，模擬區降雨量多年平均值為605 mm，水面蒸發量多年平均值為592.9 mm。由于模擬區位于渭河一級階地，根據西安市地貌圖及模擬區相關地質資料[21-22]，確定模擬區降雨入滲補給系數為0.37，在模型中通過PRECHARGE模塊進行設置。(2)對于地下水側向徑流補給，運用MODFLOW的Wells模塊，在模型中設置抽注水井，通過側向補給量確定井的流量。

4.3.4 模擬時長及初始流場 模擬期為2017年3月-2022年3月，在設置好模型各項參數后進行地下水穩定流場的模擬，得到模型的初始流場，如圖8所示。

4.4 模型識別和校正

選擇地下水流數值模擬中常用的“試錯法”[19]對模型參數進行校正，即不斷地對比觀測井的模擬計算水位和實測水位，根據對比結果針對性地調整水文地質參數，直至觀測井的模擬計算水位和實測水位誤差在合理范圍之內，則認為模型可用于區域流場的預測。在建立確定區域水文地質參數分區的過程中，參考斗門水庫建設可行性研究結果，充分收集斗門水庫試驗段建設前期工程勘測資料[22-23]，劃分出模型的水文地質參數分區。在此基礎上對模型進行水文地質參數校正：選取2017年3月-2018年3月為模型的識別和校正期，根據2018年3月庫岸邊監測井的實測地下水位繪制校正對比流場，從而得出校正后的各分區及其對應的參數值，如圖9、表3所示。

表3 模擬區各分區水文地質參數表

圖5 模擬區范圍示意圖 圖6 模擬區網格剖分示意圖

圖7 模擬區頂、底板高程(單位：m)

圖8 模擬區初始流場(單位：m) 圖9 模擬區水文地質參數分區圖

4.5 結果與分析

利用模型識別和校正后的水文地質參數和邊界條件，采用有限差分軟件Visual MODflow 模擬了斗門水庫蓄水5年后的地下水流場變化，模擬結果如圖10所示。

圖10 斗門水庫蓄水5年后地下水流場(單位：m)

通過蓄水前后(圖8及圖10)的地下水流場模擬結果對比可以看出：蓄水前，受地形因素控制，模擬區地下水流向為東南向西北，整體地下水位介于387～396 m之間，庫區地下水位介于391～394 m之間；蓄水5年后，模擬區地下水位抬升約0.1～0.5 m，整體地下水位介于388.5～396 m之間，其中庫區東北部水位抬升較為明顯，等水位線變得更為密集，水力梯度變大，地下水流向西發生偏轉，結合表1、表2分析其原因是由于庫區東北部滲漏強度較大，對原有的地下水流場產生了較大影響；庫區南部地下水流場變化較小，地下水流向基本保持不變。綜合以上分析可知，水庫蓄水會對區域地下水流場產生一定影響，使得區域地下水位尤其是庫區東北部地下水位產生了較為明顯的抬升，地下水流向發生偏轉，應當加強庫區東北部的的地下水位動態監測，做好防滲措施。

5 結 論

(1)斗門水庫試驗段庫底沉積物巖性以淤泥質土和壤土為主，東北部地層中夾有少量中砂；庫底各分區滲透系數差異較小，量級均為10-3m/d，庫區平均滲漏強度約為0.041 mm/d，即14.97 mm/a；庫區東北部沉積物滲透系數及滲漏強度均較大，應當加強庫區東北側的地下水位動態監測，做好庫底防滲措施。

(2)庫區平均滲漏量約為19.37 m3/d，即7070.05 m3/a，年滲漏損失約占試驗段總庫容的4.6‰；庫區滲漏量受降水量影響較為明顯，豐水季節的滲漏量較小。

(3)水庫蓄水對周邊區域原有的地下水流場產生了一定影響，蓄水5年后使得模擬區地下水位整體抬升0.1～0.5 m左右，庫區東北部水力梯度變大，地下水流向發生偏轉。