李文溢， 楊阿敏， 周維博

(1.長安大學(xué) 環(huán)境科學(xué)與工程學(xué)院， 陜西 西安 710054；2.長安大學(xué) 旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點實驗室， 陜西 西安 710054)

1 研究背景

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

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

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

本文以陜西省西咸新區(qū)斗門水庫試驗段為例，采用數(shù)值法和解析法相結(jié)合，分析庫區(qū)滲漏量及水庫蓄水對區(qū)域地下水流場的影響，以期為庫底防滲、水庫調(diào)蓄、地下水環(huán)境保護及斗門水庫的進一步建設(shè)提供相關(guān)科學(xué)依據(jù)。

2 研究區(qū)概況

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

圖1 研究區(qū)地理位置圖

2.1 地形地貌

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

2.2 地層巖性

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

2.3 水文地質(zhì)條件

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

3 庫區(qū)滲漏量計算

3.1 水庫滲漏模型

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

(1)

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

圖2 水庫滲漏模型示意圖

3.2 水庫滲漏分區(qū)

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

圖3 庫區(qū)滲漏分區(qū)圖

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

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

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

(2)

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

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

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

表1 滲漏分區(qū)相關(guān)參數(shù)

3.3 水庫滲漏計算

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

表2 庫區(qū)滲漏量及滲漏強度

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

4 地下水流數(shù)值模擬

4.1 水文地質(zhì)概念模型的建立

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

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

4.2 地下水流系統(tǒng)數(shù)學(xué)模型

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

(3)

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

4.3 模型設(shè)計

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

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

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

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

4.4 模型識別和校正

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

表3 模擬區(qū)各分區(qū)水文地質(zhì)參數(shù)表

圖5 模擬區(qū)范圍示意圖 圖6 模擬區(qū)網(wǎng)格剖分示意圖

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

圖8 模擬區(qū)初始流場(單位：m) 圖9 模擬區(qū)水文地質(zhì)參數(shù)分區(qū)圖

4.5 結(jié)果與分析

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

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

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

5 結(jié) 論

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

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

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