平原水庫滲流場模擬及滲漏量計(jì)算
——以斗門水庫試驗(yàn)段為例

李文溢， 楊阿敏， 周維博

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

1 研究背景

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

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

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

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

2 研究區(qū)概況

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

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

2.1 地形地貌

斗門水庫試驗(yàn)段所處地貌單元為秦嶺北麓近山前沖洪積平原，庫區(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，受到人類活動(dòng)影響，頂部人工堆積零星分布(約0.5 m，主要為水泥塊、磚塊等建筑垃圾)，于庫區(qū)東北部(渭河一級階地)分布。

2.3 水文地質(zhì)條件

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

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

3.1 水庫滲漏模型

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

(1)

式中：Q為庫區(qū)總滲漏量，即庫水補(bǔ)給地下水的總量,m3/d；Qi為分區(qū)滲漏量，即庫水在各個(gè)分區(qū)對地下水的補(bǔ)給量,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ū)個(gè)數(shù)。

圖2 水庫滲漏模型示意圖

3.2 水庫滲漏分區(qū)

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

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

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

圖4 變水頭滲透試驗(yàn)法示意圖

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

(2)

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

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

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

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

3.3 水庫滲漏計(jì)算

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

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

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

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

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

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

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

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

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

(3)

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

4.3 模型設(shè)計(jì)

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

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

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

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

4.4 模型識別和校正

選擇地下水流數(shù)值模擬中常用的“試錯(cuò)法”[19]對模型參數(shù)進(jìn)行校正，即不斷地對比觀測井的模擬計(jì)算水位和實(shí)測水位，根據(jù)對比結(jié)果針對性地調(diào)整水文地質(zhì)參數(shù)，直至觀測井的模擬計(jì)算水位和實(shí)測水位誤差在合理范圍之內(nèi)，則認(rèn)為模型可用于區(qū)域流場的預(yù)測。在建立確定區(qū)域水文地質(zhì)參數(shù)分區(qū)的過程中，參考斗門水庫建設(shè)可行性研究結(jié)果，充分收集斗門水庫試驗(yàn)段建設(shè)前期工程勘測資料[22-23]，劃分出模型的水文地質(zhì)參數(shù)分區(qū)。在此基礎(chǔ)上對模型進(jìn)行水文地質(zhì)參數(shù)校正：選取2017年3月-2018年3月為模型的識別和校正期，根據(jù)2018年3月庫岸邊監(jiān)測井的實(shí)測地下水位繪制校正對比流場，從而得出校正后的各分區(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ū)地下水流向?yàn)闁|南向西北，整體地下水位介于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ū)東北部滲漏強(qiáng)度較大，對原有的地下水流場產(chǎn)生了較大影響；庫區(qū)南部地下水流場變化較小，地下水流向基本保持不變。綜合以上分析可知，水庫蓄水會(huì)對區(qū)域地下水流場產(chǎn)生一定影響，使得區(qū)域地下水位尤其是庫區(qū)東北部地下水位產(chǎn)生了較為明顯的抬升，地下水流向發(fā)生偏轉(zhuǎn)，應(yīng)當(dāng)加強(qiáng)庫區(qū)東北部的的地下水位動(dòng)態(tài)監(jiān)測，做好防滲措施。

5 結(jié) 論

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

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

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