基于Modflow和LAK3的三泉水庫滲漏量研究

王亞琴 楊軍耀 劉慶 任柳妹

摘要：水庫滲漏補給是泉域復(fù)流的重要形式，確定水庫滲漏量及影響因素對于預(yù)測泉域復(fù)流很有必要。以三泉水庫為例，運用地下水模擬軟件CMS中的Modflow和LAK3程序包，建立水庫和地下水含水層耦合模型，模擬預(yù)測不同關(guān)井壓采方案對水庫滲漏量及古堆泉口水位的影響程度，并結(jié)合正交試驗選取滲漏區(qū)域的垂向滲透系數(shù)、平均滲透性能、給水度和貯水系數(shù)進行參數(shù)敏感性分析。結(jié)果表明：實施關(guān)井壓采后，水庫滲漏量增加，泉口水位恢復(fù)加快，對滲漏量的影響大小順序為平均滲透性能>給水度>垂向滲透系數(shù)>貯水系數(shù)。

關(guān)鍵詞：水庫滲漏量;GMS;參數(shù)敏感性;三泉水庫;古堆泉

中圖分類號：TV697.2; P641.2

文獻標志碼：A

doi：10.3969/j.issn.1000- 1379.2019.02.013

1 研究區(qū)概況

古堆泉域地處山西省南部，地理坐標為東經(jīng)111°5′-112°5′，北緯35°28′-36°0 ′。古堆泉出露于整個泉域的西側(cè)，位于運城市新絳縣古堆村，主要含水層為寒武、奧陶系灰?guī)r及其上覆第四系松散沉積物層。受區(qū)域構(gòu)造影響，古堆泉成為山西省重要的巖溶大泉之一。隨著社會經(jīng)濟的發(fā)展，受井群開采和礦坑排水等影響，古堆泉水出現(xiàn)干枯斷流。

三泉水庫也位于泉域的西側(cè)，壩址臨近古堆泉。該水庫在試通水期間，水位下降，水庫滲漏量較大，庫尾滲漏區(qū)出現(xiàn)了18處大小不等的滲水塌陷坑，塌陷坑和主要泉眼高程為442.0～ 443.6 m。經(jīng)滲漏分析，在覆蓋層較厚的區(qū)域，庫水沿著原有的泉水出溢通道回灌補給巖溶水。主滲漏區(qū)位于水庫庫尾，即古堆泉口，水庫滲漏量的大小直接影響著泉域的巖溶地下水量及古堆泉口水位變化。

水庫滲漏作為泉域復(fù)流的一種方式，研究水庫滲漏量對于預(yù)測泉域復(fù)流具有重要的現(xiàn)實意義。水庫滲漏量的計算方法有很多，較為成熟方便的是利用軟件模擬計算滲漏量[1-4]。筆者運用地下水模擬軟件GMS中的Modflow和IA K3程序包，建立庫水和地下水耦合模型，預(yù)測水庫在不同運行水位并施行關(guān)井壓采情況下的滲漏量變化及對泉域巖溶水位的影響，通過監(jiān)測泉口水位計算古堆泉復(fù)流年限，并結(jié)合正交試驗來分析水庫滲漏量的主要影響因素。

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

目前應(yīng)用廣泛的地下水數(shù)值模擬軟件GMS不僅擁有良好的可視化操作界面，而且運算效率高，處理功能強大[5]。本研究選用GMS中的Modflow和LAK3程序包進行三泉水庫滲漏過程模擬。首先根據(jù)泉域的水文地質(zhì)條件建立水文地質(zhì)概念模型，然后運用Modflow程序建立地下水流模型，將水庫條件通過LAK3模塊添加進去，最后將庫水與地下水耦合起來計算。

2.1 水文地質(zhì)概念模型

（1）計算區(qū)域及邊界概化。模擬區(qū)地下水由北東、東南方向向古堆泉方向流動，為保證東部邊界水位變化不受三泉水庫滲漏的影響，設(shè)定475 m巖溶地下水等水位線作為東部邊界，概化為定水頭邊界：北部以汾陽嶺北緣斷裂為邊界，該斷層以北巖溶含水層埋深在800 m以上，巖溶地下水處于循環(huán)緩慢的滯留性狀態(tài)，故將該邊界概化為滯留性隔水邊界：西部九原山一汾陽嶺西緣斷裂概化為隔水邊界：南部侯馬盆地南緣斷裂概化為隔水邊界;模擬區(qū)計算面積約1 026 km。

（2）含水層概化。第一層為上部廣泛分布的第四系孔隙水含水層和局部分布的巖溶水、風化裂隙水含水層，底板埋深30～350 m;第二層為第三系弱透水層和巖漿巖極弱透水層，底板埋深為175—495 m;第三層為下部的奧陶系巖溶水含水層。

（3）源匯項。研究區(qū)內(nèi)的源匯項主要由5部分組成：①降水人滲補給;②灌溉回滲補給（主要為引黃灌溉水）;③三泉水庫滲漏補給（滲漏主要發(fā)生于水庫庫尾）;④蒸發(fā)（只考慮汾澮河谷范圍蒸發(fā)）;⑤人工開采（模擬區(qū)內(nèi)的農(nóng)業(yè)灌溉水井和農(nóng)村生活用水井均概化為面井，開采量采用估值并在模型識別時調(diào)整）。

2.2 數(shù)值模型

（1）數(shù)學(xué)模型建立。研究區(qū)內(nèi)含水層水文地質(zhì)參數(shù)不均一，位于泉口附近的水庫滲漏明顯，故按非均質(zhì)、各向異性三維非穩(wěn)定流進行數(shù)值模擬，數(shù)學(xué)模型為

（2）水庫（LAK3）模塊。水庫通過LAK3模塊概化，利用Triangle Interplation方法進行插值生成水庫TIN，從而確定水庫的形狀和位置。需要設(shè)置的參數(shù)主要有初始水位、水庫平均滲透性能、降水量、徑流量、供水量，其中：水庫平均滲透性能K’由水庫庫底淤積物滲透系數(shù)和淤積物厚度比值決定，水庫水位根據(jù)水庫實際運行方式賦值。

（3）參數(shù)分區(qū)。根據(jù)水文地質(zhì)特征、地形地貌單元及地層巖性將模擬區(qū)劃分為汾澮河谷階地區(qū)（I）、碳酸鹽巖深埋區(qū)（Ⅱ）、碳酸鹽巖淺埋區(qū)（Ⅲ）、碳酸鹽巖裸露區(qū)（Ⅳ）、塔兒山一九原山中低山區(qū)cv）s個參數(shù)分區(qū)，其中：中低山區(qū)根據(jù)巖性差異性又可分為碳酸鹽巖一巖漿巖淺覆蓋區(qū)（ vi）和碳酸鹽巖裸露區(qū)（v2）兩個亞區(qū)（見圖1）。

（4）模型求解。根據(jù)模擬區(qū)水文地質(zhì)條件、三泉水庫位置和模擬區(qū)巖溶含水層埋深狀況，將模型剖分為109x105x3（列×行×層）共18 138個計算單元（活動單元）。模擬期為2016年6月-2056年1月，模擬總時長為39.5 a。

（5）模型校正和驗證。選取2016年6-9月作為校正期，對模型中的參數(shù)進行校正。選取2016年10月-2017年1月作為模型驗證期，運用該階段觀測孔資料來驗證模型。經(jīng)過計算水位和實測水位的多次擬合，觀測孔的擬合誤差均小于10%，符合《地下水資源管理模型工作要求》[6]對模型誤差的要求，說明該模型參數(shù)取值合理。校正末期巖溶水位等值線見圖2。

3 滲漏量與水位預(yù)測

3.1 方案預(yù)設(shè)

為研究水庫不同蓄水位狀態(tài)，結(jié)合不同關(guān)井壓采方案的水庫滲漏量以及古堆泉口巖溶水位變化特征，設(shè)定4種方案：①現(xiàn)狀開采條件下三泉水庫低水位運行，運行水位為443.0 m;②現(xiàn)狀開采條件下水庫正常調(diào)節(jié)庫水位運行，運行水位為447.5 m，③水庫正常蓄水運行下實施一期關(guān)井壓采方案：④水庫正常蓄水運行下同步實施一二期關(guān)井壓采方案。其中：一期關(guān)井壓采方案為2017-2018年壓采巖溶水342.9萬m：二期關(guān)井壓采方案為2020-2021年壓采巖溶水706.9萬m。

3.2 預(yù)測分析

（1）模擬期內(nèi)，水庫滲漏影響最為顯著的古堆泉口水位變化在4種方案下均經(jīng)歷3個階段，即短期快速抬升一多起伏波動抬升一等幅抬升（見圖3）。古堆泉口水位由379.7 m恢復(fù)至443.0 m（即復(fù)流），4種方案下分別需要39、39、37、35 a。

（2）水庫低水位和正常蓄水位條件下，末期滲漏量相差不大，均約為1.105 0萬m/d，初期的水庫滲漏量在正常蓄水位下為2.183 8萬m/d，較低水位下大0.368 2萬m/d。水庫正常蓄水位條件下實施一期關(guān)井壓采和同步實施一二期關(guān)井壓采方案后的末期滲漏量相差不大，較正常蓄水位條件下分別大0.002 0萬、0.002 1萬m/d。4種方案下滲漏量變化趨勢相似，以水庫低水位運行條件下為例，水庫滲漏量變化和泉口水位變化歷時曲線見圖4，隨泉口水位上升水庫滲漏量迅速變小，隨后波動上升并逐漸趨于穩(wěn)定。

（3）水庫蓄水之后，庫水通過泉口滲漏通道人滲，滲漏速度快，滲漏量大;地下水漏斗具有中間小、四周大的形態(tài)特征，古堆泉水位在開始時恢復(fù)速度相對較快;隨著水庫淤積加深，淤積厚度增大，滲漏速度變慢，滲漏量減少：庫水滲漏補給地下水漏斗，隨著泉口水位抬高，相應(yīng)的水頭壓力變大，漏斗面積增大，水位回升速度相對變慢：隨著水位持續(xù)回升，在水庫滲漏補給古堆泉的過程中，含水層結(jié)構(gòu)發(fā)生變化，包氣帶變?yōu)轱査畮?，水庫相當于頂托式地往漏斗?nèi)灌水，滲漏量相對增大，滲漏速度也變快;最后水位恢復(fù)到泉口時，庫水與巖溶水達到相同水頭，滲漏量和水位回升達到平衡。

4 參數(shù)敏感性分析

4.1 介質(zhì)參數(shù)及靈敏度指標選取

在地下水流模型建立過程中，選取與巖溶水流場相關(guān)性較強的因素作為敏感性分析因子，模擬在水庫正常蓄水條件下同步實施一二期關(guān)井壓采方案的水庫滲漏情況。選定水庫所在分區(qū)的垂向滲透系數(shù)、水庫平均滲透性能、給水度、貯水系數(shù)4個因素作為敏感性分析因子，同時選出5組水平值得到四因子五水平的正交設(shè)計表L25（5）[7]。在設(shè)計正交試驗表時，將以上選取因素作同等幅度的變化，以便在相同條件下比對參數(shù)的敏感度。選取模型最終識別值為標準值，其余水平值在該值基礎(chǔ)上分別上下調(diào)整20%、40%。選取水庫的初期滲漏量作為評價指標。各參數(shù)水平值見表1。

4.2 試驗結(jié)果分析

將這25組試驗值代人模型中進行GMS計算，利用在地下水數(shù)值模擬研究中創(chuàng)建專用圖形的GW_Chart程序的運算器和水文提取器功能輸出LakePolts.從而得到滲漏量隨時間的變化值。正交試驗結(jié)果及初期滲漏量見表2。

最常用的正交試驗數(shù)理統(tǒng)計分析方法是極差分析法，即通過極差的大小來反映評價指標對各個因素的敏感程度，極差越大說明指標對該因素的敏感性越高，反之則敏感性越低[8]。對正交試驗結(jié)果采用極差法進行分析，結(jié)果見表3。影響水庫滲漏量的因素極差大小順序為平均滲透性能>給水度>垂向滲透系數(shù)>貯水系數(shù)。

通過計算4個因子在不同水平下的極差可知，水庫初期滲漏量對平均滲透性能最為敏感，其次為給水度，而對滲漏區(qū)域的垂向滲透系數(shù)和貯水系數(shù)敏感性較低，因此在數(shù)值模擬過程中對水庫的平均滲透性能準確性要求最高。

為研究水庫平均滲透性能K’與水庫初期滲漏量Q的關(guān)系，在其他條件不變的情況下只改變K’的大小，得到K’與水庫初期滲漏量Q的關(guān)系曲線（見圖5）?？芍?，水庫平均滲透性能K’與初期滲漏量Q基本成線性正相關(guān)關(guān)系，表達式近似為Q=1 030.7K'+ 0.086 7，即水庫平均滲透性能K’越大，初期滲漏量Q也越大。

5 結(jié)論

（1）數(shù)值模擬結(jié)果表明，水庫低水位和正常蓄水位條件下，水庫末期滲漏量相差不大，均約為1.105 0萬m/d;實施關(guān)井壓采后，方案③和方案④末期滲漏量較正常蓄水位條件下分別大0.002 0萬、0.0021萬m/d。

（2）在方案③和方案④下，泉口水位恢復(fù)至443 m（即復(fù)流）分別需要37、35 a，較未實施關(guān)井壓采方案分別提前2、4a，說明關(guān)井壓采力度加大后可提高泉口水位的恢復(fù)速度。

（3）參數(shù)敏感性分析表明，影響水庫滲漏量各因素極差大小順序為平均滲透性能>給水度>垂向滲透系數(shù)>貯水系數(shù)，其中：水庫的平均滲透性能K’值與初期滲漏量Q成線性正相關(guān)關(guān)系，表達式近似為Q=1030.7K’+0.086 7。因此在地下水流隨機模擬過程中，應(yīng)對水庫平均滲透性能的準確性提出更高要求。

參考文獻：

[1]程春龍，束龍倉，魯程鵬，等，應(yīng)用MODFLOW和LAK3計算水庫滲漏量[J].工程勘察，2013，41（7）：31-34.

[2] 呂路，楊軍耀，秦嘉楠，汾河三期蓄水工程人工湖滲漏對太原市淺層地下水影響的數(shù)值模擬[J].水電能源科學(xué)，2017，35（1）：89 - 92.

[3]介飛龍，朱志新，封麗華，等.LAKE在水庫地下水影響預(yù)測中的應(yīng)用[J].中國水運，2016，16（1）：121-123.

[4]劉涵宇，夏強，張世殊，等，基于水文地質(zhì)參數(shù)敏感度分析的水庫滲漏量不確定性評價[J].科學(xué)技術(shù)與工程，2017，17（14）：32-38.

[5]祝曉彬，地下水模擬系統(tǒng)（GMS）軟件[J].水文地質(zhì)工程地質(zhì)，2003（5）：53-55.

[6] 王兆馨，林學(xué)鈺，余國光，等，地下水資源管理模型工作要求：GB/T 14497-1993[S].北京：中國標準出版社，1994：8.

[7]劉瑞江，張業(yè)旺，聞崇煒，等，正交試驗設(shè)計和分析方法研究[J].實驗技術(shù)與管理，2010，27（9）：52-55.

[8] 郝靜，賈仰文，張永祥，等，應(yīng)用正交試驗法分析地下水流模型參數(shù)靈敏度[J].人民黃河，2015，37（9）：66-68.