李鳳麗，徐嘉璐，張 游，韓 宇

(1.水發(fā)規(guī)劃設(shè)計(jì)有限公司，濟(jì)南 250013；2.山東省湖泊流域管理信息化工程技術(shù)研究中心，濟(jì)南 250013)

0 引 言

地下水庫(kù)起源于地下水人工補(bǔ)給，可以通過(guò)天然地下儲(chǔ)水空間攔蓄、調(diào)節(jié)和利用地下水[1, 2]。1972年，日本松尾氏采用灌漿法在長(zhǎng)崎縣野母崎町樺島上建設(shè)了世界上第一座地下水庫(kù)，隨后又在福井、福岡、沖繩縣等地建設(shè)了十余座地下水庫(kù)，不僅提高了工、農(nóng)業(yè)供水能力，還可以防止海水入侵[3, 4]。20世紀(jì)80年代以來(lái)，美國(guó)、瑞典、荷蘭和德國(guó)等地逐漸開(kāi)展了人工儲(chǔ)水與回采工程，先將地表水回灌到地下水庫(kù)中，再通過(guò)排水渠和抽水井提取地下水，對(duì)水質(zhì)進(jìn)行適當(dāng)處理后供城市用水[5, 6]。我國(guó)第一座地下水庫(kù)于1975年在河北省南宮市興建[7]，之后北京西郊、山東龍口、大連旅順、廣西、貴州和福建等地[8-12]先后建設(shè)了地下水庫(kù)，為利用地下水提供了寶貴的經(jīng)驗(yàn)。在地下水庫(kù)水流數(shù)值模擬方面[13,14]，藍(lán)盈盈[15]運(yùn)用GMS軟件對(duì)洮兒河扇形地地下水庫(kù)水流進(jìn)行了數(shù)值模擬，結(jié)果表明模型可以用于地下水位的預(yù)報(bào)；中科院利用回灌試驗(yàn)和數(shù)值模擬研究了北京市西郊地下水庫(kù)永定河河道的入滲能力[16]；為了研究秦嶺山前洪積扇地下水庫(kù)調(diào)蓄功能，康華等[17]在太平河洪積扇進(jìn)行回灌試驗(yàn)，并使用Visual Modflow進(jìn)行數(shù)值模擬，結(jié)果表明適度加強(qiáng)開(kāi)采可以激發(fā)地下水庫(kù)更大的調(diào)蓄潛力。

威海市區(qū)和文登區(qū)屬于資源性嚴(yán)重缺水的地區(qū)，隨著社會(huì)經(jīng)濟(jì)的發(fā)展，水資源供需矛盾日益加劇，為了緩解這一矛盾，在母豬河下游擬建地下水庫(kù)，進(jìn)行河道攔蓄，與位于其上游的米山水庫(kù)進(jìn)行聯(lián)合調(diào)度，增加地表水和地下水的水資源利用率。本文結(jié)合母豬河下游的水文地質(zhì)條件，利用Visual Modflow進(jìn)行數(shù)值模擬，對(duì)地下水庫(kù)進(jìn)行調(diào)蓄分析并選定出較為合適的開(kāi)采方案。

1 研究區(qū)概況

威海市位于山東半島東端，地處北緯36°41′～ 37°35′、東經(jīng)121°11′～122°42′。北、東、南三面瀕臨黃海，西與煙臺(tái)市接壤。威海市區(qū)、文登區(qū)水資源有年際變幅大、年內(nèi)分配不均、地域分布不均、開(kāi)發(fā)難度大、利用率低、人均水資源量少等特點(diǎn)，人均水資源占有量為858 m3，僅為全國(guó)人均占有量的32%，屬于資源性嚴(yán)重缺水地區(qū)。近年來(lái)隨著經(jīng)濟(jì)社會(huì)的發(fā)展，各部門用水量急劇增加，水資源供需矛盾日益嚴(yán)重。研究開(kāi)發(fā)新的供水工程，可增加威海市區(qū)、文登區(qū)供水能力，緩解水資源短缺狀況。

母豬河是威海市文登區(qū)第一大河，干流總長(zhǎng)65 km，流域面積1 115.18 km2，占文登區(qū)總面積的63%，流域圖見(jiàn)圖1。母豬河中上游建有米山水庫(kù)，是一座以防洪、灌溉、城市及工業(yè)供水為主，兼顧發(fā)電和養(yǎng)殖等綜合效益的多年調(diào)節(jié)大(2)型水庫(kù)，控制流域面積440 km2，占母豬河流域面積的39.46%。母豬河地下水水庫(kù)擬建在母豬河下游，下游河段河谷寬闊，下游河床、河漫灘上伏6.0～10.0 m厚的第四系覆蓋層，以壤土及中粗砂為主構(gòu)成，局部有中細(xì)砂、細(xì)砂，屬中等強(qiáng)透水層，地下水埋深一般為1.0～3.5 m，下伏基巖為二長(zhǎng)花崗巖，地下截滲條件較好，水文地質(zhì)剖面圖見(jiàn)圖2。母豬河地下水庫(kù)不僅可以發(fā)揮多年調(diào)節(jié)作用，還能攔蓄地表徑流用以補(bǔ)給地下水，從而增加地下水儲(chǔ)量和開(kāi)采量。此外，文登區(qū)雖然有較多零星水系，但是大多獨(dú)立入海，無(wú)法攔蓄，使得水資源可利用率降低，且母豬河下游由于地下水長(zhǎng)期超采，形成地下水位降落漏斗，在河口誘發(fā)海水入侵。興建地下水庫(kù)工程，通過(guò)“一截一蓄”，進(jìn)行河道攔蓄，增加地表水的利用量；通過(guò)“以豐補(bǔ)歉”，在豐水期蓄積大量的水資源以備枯水期利用，還可以為米山水庫(kù)提供應(yīng)急備用水源，緩解區(qū)域水資源供需矛盾；通過(guò)“蓄淡壓咸”，防止海水入侵，改善區(qū)域生態(tài)環(huán)境。

圖1 母豬河流域圖(單位：m)Fig.1 The basin map of Muzhu River

圖2 母豬河水文地質(zhì)剖面圖Fig.2 The hydrogeological profile of Muzhu River

文登區(qū)母豬河地下水庫(kù)工程分為上級(jí)庫(kù)工程與下級(jí)庫(kù)工程，主要建筑物包括地下截滲墻、橡膠壩等。地下水庫(kù)庫(kù)容和水位參數(shù)見(jiàn)表1。上級(jí)地下截滲工程位于院東村東側(cè)，共長(zhǎng)1.9 km，頂高程為9.0 m。下級(jí)地下截滲工程位于新建高速路南側(cè)，長(zhǎng)共1.8 km，頂高程為6.0 m。該工程擬攔蓄母豬河右岸一級(jí)支流河水及兩岸的地下水。母豬河上下級(jí)庫(kù)示意圖見(jiàn)圖3。

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

2.1 確定模擬范圍

考慮到降雨蒸發(fā)的時(shí)空分布對(duì)庫(kù)區(qū)的影響，為使模擬結(jié)果更加符合實(shí)際情況，以分水嶺為邊界對(duì)庫(kù)區(qū)的整體匯流區(qū)域進(jìn)行模擬。圖4中紅線為母豬河地下水庫(kù)流域分水嶺，黃線為米山水庫(kù)流域分水嶺。由于區(qū)域②的匯流被米山水庫(kù)所截留，因此母豬河地下水庫(kù)匯流面積僅為①區(qū)域，選?、賲^(qū)域作為模擬范圍。

表1 母豬河地下水庫(kù)參數(shù)Tab.1 Parameters of underground reservoir of Muzhu River

圖3 母豬河上下級(jí)庫(kù)示意圖Fig.3 Schematic diagram of the upper and lower reservoirs of Muzhu River

圖4 流域水系圖Fig.4 The drainage map of basins

2.2 含水層結(jié)構(gòu)概化

根據(jù)地下水埋藏條件，母豬河地下水庫(kù)模擬區(qū)含水層結(jié)構(gòu)確定為潛水含水層。綜合考慮各地層的水力聯(lián)系程度、含水層滲透性大小以及巖性特征，對(duì)含水層結(jié)構(gòu)自上而下進(jìn)行概化。研究區(qū)地表以及下伏土壤多為中粗砂且滲透性相近，本次模擬將其概化為一層，該層滲透性強(qiáng)，富水性好，為母豬河地下水庫(kù)主要富水層，底層為強(qiáng)風(fēng)化層花崗巖，透水率低，概化為不透水層。

2.3 邊界條件概化

通過(guò)對(duì)母豬河地下水庫(kù)地下水系統(tǒng)的分析，發(fā)現(xiàn)母豬河地下水庫(kù)為河道型地下水庫(kù)，其儲(chǔ)水地層均位于母豬河及支流河床上。其邊界條件概化如下：①底邊界：庫(kù)區(qū)基巖為下元古界荊山群變質(zhì)巖，透水性較弱，富水性差，隔水性能良好，概化為隔水層，為零通量邊界。②上邊界：上邊界為潛水面，巖性主要為富水性和滲透性強(qiáng)的中粗砂，受大氣降雨、人為開(kāi)采和蒸發(fā)排泄等因素的影響，為變通量邊界。③東部邊界和西部邊界：根據(jù)流域分水嶺的劃分東、西部邊界，均為隔水邊界。④流域上游邊界：上游邊界米山水庫(kù)以東為分水嶺設(shè)置為隔水邊界，米山水庫(kù)根據(jù)實(shí)際放流情況和地下水測(cè)流補(bǔ)給程度，設(shè)置為變通量邊界條件。⑤流域下游邊界：下游地勢(shì)較低為母豬河地下水庫(kù)的主要排泄通道，根據(jù)水位和地層高程，設(shè)置為自由排水邊界條件。

2.4 水文地質(zhì)參數(shù)獲取和分區(qū)

用于地下水流模型的水文地質(zhì)參數(shù)主要分為兩類：第一類是用于計(jì)算地下水源匯項(xiàng)的經(jīng)驗(yàn)參數(shù)，如大氣降水入滲系數(shù)、灌溉入滲系數(shù)、蒸發(fā)系數(shù)等；第二類是含水層的水文地質(zhì)參數(shù)，主要包括潛水含水層的滲透系數(shù)(K)、給水度(Sy)和有效孔隙度(η)等。在建模工作中，首先需要根據(jù)水文地質(zhì)條件等因素確定參數(shù)的分區(qū)。經(jīng)分析發(fā)現(xiàn)母豬河地下水庫(kù)流域的地貌特征主要分為丘陵和平原兩種，丘陵區(qū)巖性主要為下元古界荊山群變質(zhì)巖，平原區(qū)巖性主要為中粗砂，二者參數(shù)差別較大。因此，模型分為兩個(gè)參數(shù)區(qū)域。圖5為母豬河地下水庫(kù)流域活動(dòng)單元的三維圖，地勢(shì)較高的深色區(qū)域?yàn)榍鹆陞^(qū)，地勢(shì)低平的淺色區(qū)域?yàn)槠皆瓍^(qū)。根據(jù)實(shí)測(cè)鉆孔探測(cè)的底板高程和設(shè)置截滲墻的影響，確定母豬河地下水庫(kù)最大回水面積為圖中黑色線所圍區(qū)域，橘色線段為截滲墻所在位置。經(jīng)過(guò)計(jì)算模擬范圍丘陵區(qū)的面積為138.2 km2，平原區(qū)的面積為80.7 km2。

按照參數(shù)的分區(qū)給定參數(shù)初值，通過(guò)參數(shù)率定得出各參數(shù)的分區(qū)值。參數(shù)分區(qū)初始值主要通過(guò)野外抽水試驗(yàn)結(jié)合地區(qū)經(jīng)驗(yàn)值，經(jīng)過(guò)綜合分析確定，具體見(jiàn)表2。

2.5 源匯項(xiàng)確定

地下水系統(tǒng)的源匯項(xiàng)主要包括補(bǔ)給項(xiàng)與排泄項(xiàng)。補(bǔ)給項(xiàng)主要包括降雨入滲補(bǔ)給、河道側(cè)滲補(bǔ)給和山前側(cè)滲補(bǔ)給。排泄項(xiàng)主要包括地下水開(kāi)采、潛水蒸發(fā)和下游排泄。

表2 主要分區(qū)參數(shù)初始值Tab.2 The initial parameters of areas

3 地下水流數(shù)學(xué)模型

3.1 模型建立

根據(jù)模擬區(qū)地下水系統(tǒng)水文地質(zhì)概念模型，建立了地下水流數(shù)學(xué)模型。當(dāng)不考慮弱透水層貯水能力時(shí)，淺層潛水系統(tǒng)地下水運(yùn)動(dòng)的數(shù)學(xué)模型概化為三維不穩(wěn)定流動(dòng)系統(tǒng)，其數(shù)學(xué)模型如下：

式中：K為邊界法線方向滲透系數(shù)，m/d；Sy為給水度；S1、S2為模擬區(qū)域第一類和第二類邊界；H0為含水層初始水頭，m；q為含水層第二類邊界單位面積過(guò)水?dāng)嗝嫜a(bǔ)給流量，m3/d；φ為源匯項(xiàng)強(qiáng)度(包括開(kāi)采強(qiáng)度、入滲強(qiáng)度、蒸發(fā)強(qiáng)度等)，1/d；D為滲流區(qū)域；n為滲流區(qū)邊界的單位外法線方向，矢量；W為源匯項(xiàng)，流進(jìn)為正，流出為負(fù)；t為時(shí)間。

3.2 網(wǎng)格剖分

此次模擬采用矩形網(wǎng)格對(duì)滲流區(qū)進(jìn)行離散化(剖分)，將復(fù)雜的滲流問(wèn)題處理成符合剖分單元內(nèi)規(guī)則的簡(jiǎn)單滲流問(wèn)題。網(wǎng)格剖面見(jiàn)圖6。模型共剖分單元格100×100(行×列)，單個(gè)單元格的尺寸是280 m×254 m，由于模擬區(qū)含水層結(jié)構(gòu)為潛水含水層，此次網(wǎng)格剖分在垂直方向上分為1層。

圖6 網(wǎng)格剖分圖Fig.6 The grid sectional map

3.3 參數(shù)離散化

根據(jù)收集的研究區(qū)地下水系統(tǒng)資料，本次數(shù)值模型的模擬期為1979年1月到1979年12月，將整個(gè)模擬期劃分為12個(gè)應(yīng)力期，每個(gè)應(yīng)力期為一個(gè)自然月，時(shí)間步長(zhǎng)為十天。在每個(gè)應(yīng)力期中，所有外部源匯項(xiàng)的強(qiáng)度保持不變。初始條件：根據(jù)前期觀測(cè)孔的監(jiān)測(cè)資料，插值計(jì)算潛水含水層的滲流場(chǎng)，并將其作為研究區(qū)的初始流場(chǎng)。模型潛水含水層初始流場(chǎng)見(jiàn)圖7。邊界條件：將所計(jì)算出來(lái)的各邊界流入流出量輸入到模型之中，通過(guò)邊界附近流場(chǎng)的擬合程度，適當(dāng)調(diào)整邊界流入流出量。

圖7 初始流場(chǎng)圖Fig.7 The initial flow field map

3.4 參數(shù)率定與驗(yàn)證

選取葛家鎮(zhèn)小英村東21號(hào)井、澤頭鎮(zhèn)獸醫(yī)站院內(nèi)25號(hào)井、葛家鎮(zhèn)糧食所院內(nèi)26號(hào)井、葛家鎮(zhèn)鋪集固恒機(jī)械廠招待所院27號(hào)井和澤頭鎮(zhèn)道口村委院內(nèi)東南角48號(hào)井共五眼井進(jìn)行參數(shù)率定，位置分布見(jiàn)圖3。采用2016年逐日地下水位資料實(shí)測(cè)值與模型模擬值進(jìn)行參數(shù)率定，結(jié)果表明實(shí)測(cè)值與模擬值擬合度較好，相關(guān)系數(shù)均大于95%。運(yùn)用pest參數(shù)迭代調(diào)整參數(shù)，調(diào)整后的平原區(qū)中粗砂滲透系數(shù)K值為38.5 m/d，而丘陵區(qū)變質(zhì)巖滲透系數(shù)K值為0.012 m/d 。平原區(qū)中粗砂的給水度為0.21，丘陵區(qū)變質(zhì)巖給水度0.058。

將調(diào)整好的參數(shù)輸入到模型中運(yùn)算，輸出結(jié)果與2017年的實(shí)測(cè)地下水位數(shù)據(jù)進(jìn)行驗(yàn)證，由于數(shù)據(jù)資料不全，只選取21號(hào)，26號(hào)和48號(hào)三眼井進(jìn)行模型驗(yàn)證，模擬值與實(shí)測(cè)值相關(guān)系數(shù)為98.7%，表明模型能夠較好的反應(yīng)實(shí)際情況，精度較高。圖8和圖9分別是模型進(jìn)行參數(shù)率定和驗(yàn)證時(shí)48號(hào)井水位擬合圖。

圖8 2016年48號(hào)井水位實(shí)測(cè)值和模擬值擬合圖Fig.8 Fitting curve of measured and simulated water level of well 48 in 2016

圖9 2017年48號(hào)井水位實(shí)測(cè)值和模擬值擬合圖Fig.9 Fitting curve of measured and simulated water level of well 48 in 2017

4 地下水庫(kù)調(diào)蓄與開(kāi)采

4.1 地下水庫(kù)調(diào)蓄分析

根據(jù)建立的模型，結(jié)合水量調(diào)算，計(jì)算出母豬河地下水庫(kù)豐水年、平水年、枯水年和特枯年的最大調(diào)蓄水量，結(jié)果見(jiàn)表3。除利用模型計(jì)算最大調(diào)蓄水量，還利用水量平衡方程對(duì)地下水庫(kù)進(jìn)行典型年調(diào)算，調(diào)算結(jié)果豐、平、枯和特枯年的分別是2 575.94、1 969.87、1 659.30和890.71 萬(wàn)m3。通過(guò)對(duì)比可知，模型模擬結(jié)果和典型年調(diào)算結(jié)果相差較小，豐水年、平水年、枯水年結(jié)果相差均在10%之內(nèi)，特枯年結(jié)果相差12.13%，再次驗(yàn)證所建模型的合理性。從水量上分析，母豬河地下水庫(kù)建成后，有較大的開(kāi)發(fā)潛力；從含水層水力條件及其富水性特征上分析，潛水含水層可調(diào)蓄空間適用性大；從長(zhǎng)遠(yuǎn)角度考慮，合理的調(diào)蓄、優(yōu)化儲(chǔ)水用水方案是合理利用地下水資源的關(guān)鍵。

4.2 地下水開(kāi)采方案選定

母豬河地下水庫(kù)建成后，擬將蓄積的地下水抽取并泵送至上游米山水庫(kù)，取代部分外調(diào)水源，優(yōu)先利用當(dāng)?shù)貎?yōu)質(zhì)的地下水資源。為合理開(kāi)采地下水，結(jié)合上游米山水庫(kù)的放水情況和母豬河地下水庫(kù)的開(kāi)采條件，初步擬定12種開(kāi)采方案，見(jiàn)表4，并通過(guò)模型分析各方案水位變化情況和疏干情況，對(duì)比每種方案的優(yōu)缺點(diǎn)。

表3 地下水庫(kù)調(diào)蓄計(jì)算表 萬(wàn)m3

表4 不同開(kāi)采方案分析結(jié)果Tab.4 The results of different scenarios of abstracting groundwater

對(duì)各方案水位變化進(jìn)行分析，從模型的穩(wěn)定性來(lái)看，流量為0.01 m3/s的方案比0.015 m3/s的方案更加穩(wěn)定；但是在開(kāi)采總量上，0.01 m3/s流量條件下的6種方案中只有55眼井全年開(kāi)采的方案接近最大調(diào)蓄水量，同等條件下0.015 m3/s流量的開(kāi)采方案能獲得較大的開(kāi)采量；從工程經(jīng)濟(jì)上分析，55眼井的方案工程量大、費(fèi)用高，而27眼井又很難滿足開(kāi)采需求。結(jié)合初步擬定的12種開(kāi)采方案的分析結(jié)果，經(jīng)過(guò)多次計(jì)算、調(diào)試和優(yōu)化，最終確定36眼開(kāi)采井(下級(jí)庫(kù)23眼、上級(jí)庫(kù)13眼)，流量為0.015 m3/s且全年開(kāi)采的方案，較為滿足工程實(shí)施和開(kāi)采量的需求。開(kāi)采井在模型中的分布如圖10所示，在開(kāi)采井的位置上布設(shè)觀測(cè)井，觀測(cè)水位的實(shí)時(shí)變化。

圖10 井的位置分布Fig.10 The location of wells

開(kāi)采井設(shè)置完畢，模擬枯水年條件下，開(kāi)采地下水對(duì)地下水位的影響情況，模擬時(shí)段末模擬區(qū)地下水降深情況如圖11所示。模型運(yùn)行過(guò)程中大多數(shù)井水位比較穩(wěn)定，表明開(kāi)采量和補(bǔ)給量達(dá)到平衡狀態(tài)，能形成穩(wěn)定的降落漏斗，在當(dāng)前開(kāi)采量下可以持續(xù)抽取地下水。但一些離河床較遠(yuǎn)的井(如5號(hào)井)，水位有較大幅度的降低，如果持續(xù)抽水可能會(huì)導(dǎo)致單元格水位疏干。因此，在實(shí)際布井的過(guò)程中，應(yīng)根據(jù)實(shí)際情況盡可能依河而建，且盡量避開(kāi)村莊和農(nóng)田。

圖11 模擬區(qū)地下水降深圖Fig.11 The drawdown of simulation zone

5 結(jié) 論

本文通過(guò)分析和計(jì)算，可以得出以下結(jié)論。

(1)模型參數(shù)率定時(shí)，實(shí)測(cè)值與模擬值相關(guān)系數(shù)均大于95%，模型驗(yàn)證時(shí)，模擬值與實(shí)測(cè)值相關(guān)系數(shù)為98.7%，表明模型能夠較好地反映實(shí)際情況。

(2)通過(guò)所建模型對(duì)地下水庫(kù)進(jìn)行調(diào)蓄分析，結(jié)果表明豐、平、枯和特枯年的最大調(diào)蓄水量分別是2 392.8、2 056.5、1 743.3和1 012.8 萬(wàn)m3。母豬河地下水庫(kù)建成后，有較大的開(kāi)發(fā)潛力，潛水含水層有較大的可調(diào)蓄空間。

(3)經(jīng)過(guò)多次計(jì)算和調(diào)試，發(fā)現(xiàn)36眼開(kāi)采井(下級(jí)庫(kù)23眼、上級(jí)庫(kù)13眼)，流量為0.015 m3/s且全年開(kāi)采的方案，能同時(shí)滿足工程要求和開(kāi)采量的需求，為較為合適的開(kāi)采方案。

本次研究成果可為威海市母豬河地下水庫(kù)建設(shè)和地下水資源開(kāi)發(fā)利用提供參考，同時(shí)，也可以為我國(guó)地下水庫(kù)建設(shè)提供借鑒和思路。