濟南市歷城區(qū)地下水數(shù)值模擬與預(yù)測

朱恒華，張曉偉，王 瑋，賈 超，劉治政，劉柏含

(1.山東省地質(zhì)調(diào)查院，山東 濟南 250013；2.中國地質(zhì)大學(xué)(武漢)，湖北 武漢 430074；3.山東大學(xué) 土建與水利學(xué)院，山東 濟南 250061； 4.山東大學(xué) 海洋研究院，山東 青島 266237)

隨著社會經(jīng)濟的不斷發(fā)展，城市工業(yè)和農(nóng)業(yè)用水量日益增加，水資源量消耗劇增。由于地表水易受污染，利用率較低，開采地下水的需求日益凸顯。然而，過度開采地下水會造成地面沉降、巖溶塌陷等地質(zhì)問題[1]。因此，科學(xué)合理地開采地下水就顯得尤為重要。近年來GMS被廣泛地應(yīng)用于地下水資源的計算評估預(yù)測及相關(guān)的環(huán)境水文地質(zhì)研究中并取得了良好的效果[2，3]。在地下水數(shù)值模型建模過程中，相關(guān)水文地質(zhì)參數(shù)對模型精度影響較大，其中又以滲透系數(shù)最為敏感[4]。PEST是GMS中自動調(diào)參模塊，因此選取滲透系數(shù)進行參數(shù)自動估計。運用PEST來調(diào)整選定的參數(shù)，重復(fù)運行MODFLOW模塊的計算，直到計算結(jié)果和野外觀測值相吻合。歷城區(qū)位于濟南市，前人研究多集中于巖溶水，對松散巖類孔隙水研究較少。本文通過GMS軟件建立歷城區(qū)淺層地下水數(shù)值模型，對不同開采方案下的地下水流場進行預(yù)報，為淺層地下水資源開采提供參考依據(jù)。

1 研究區(qū)地質(zhì)概況

歷城區(qū)位于濟南市中部，地勢南高北低，如圖1。根據(jù)區(qū)內(nèi)地貌特征，可劃分為堆積地貌和侵蝕地貌，本區(qū)南部是一個以古生代地層為主體的北傾單斜構(gòu)造，北部則是沉積了厚度較大的第四系地層和第三系松散堆積物的凹陷，區(qū)內(nèi)地層由老到新依次為古生界寒武系、奧陶系、石炭-二疊系；新生界第三系及第四系。

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

研究區(qū)處于魯中山地的北緣以及黃河下游沖洪積平原區(qū)，區(qū)內(nèi)地層、構(gòu)造條件復(fù)雜，根據(jù)地層巖性、含水介質(zhì)及地下水運動、儲存等特征，區(qū)內(nèi)地下水類型可分為松散巖類孔隙水、碎屑巖類孔隙裂隙水、碳酸鹽巖類裂隙巖溶水。除基巖裸露區(qū)外，松散巖類孔隙水在全區(qū)普遍分布。含水層巖性為粉土、粉細砂、中細砂、粗砂、礫石及粘質(zhì)砂土等，地下水多屬潛水或微承壓水，如圖2；碎屑巖類孔隙裂隙水主要分布于北部、西北部，孔隙裂隙不甚發(fā)育，富水性較差；碳酸鹽巖類裂隙巖溶水主要賦存在奧陶系灰?guī)r中，分布于南部及黃河以北局部地段，與淺層孔隙水水力聯(lián)系微弱。

圖2 黃河沖積平原區(qū)水文地質(zhì)剖面圖

2 地下水流模型

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

2.1.1 地層概化

研究區(qū)主要含水層確定為淺層孔隙水含水層，第四系厚度20～170 m，總體呈現(xiàn)由南向北漸厚的趨勢。淺層孔隙水以下為中層咸水，不進行地下水開采，故將含水層概化為深度170 m的非均質(zhì)、各向同性的潛水微承壓含水層。

2.1.2 邊界條件

1)垂向邊界：上邊界接受大氣降水、排泄方式為大氣蒸發(fā)，是不斷變化的水量交換邊界，定為流量邊界；下邊界因水量交換微弱，概化為相對隔水邊界。

2)側(cè)向邊界：研究區(qū)南部為基巖出露區(qū)，概化為隔水邊界[5]；北部以李家斷層為界概化為隔水邊界，其余邊界作為流量邊界，根據(jù)達西定律計算截面流量。

2.1.3 源匯項

(1)大氣降水補給量以面狀補給賦值，根據(jù)降雨入滲系數(shù)法[6]計算本區(qū)松散巖類孔隙水大氣降水補給總量為4 249.965萬 m3/a(合11.64萬 m3/d)。灌溉回滲量以面狀補給賦值，依據(jù)研究區(qū)內(nèi)的灌溉定額、類型及包氣帶巖性等計算區(qū)內(nèi)灌溉回滲補給總量為187.663萬 m3/a，其中井灌回滲補給量為59.251萬 m3/a，黃灌回滲補給量為132.973萬 m3/a；區(qū)內(nèi)對地下水具有一定補排規(guī)模的河渠有黃河、邢家度總干渠。河渠滲漏量以線狀補給賦值[7]。根據(jù)測流數(shù)據(jù)，黃河斷面的單側(cè)單寬補排量為1.102 m3/d·m，側(cè)滲補給總量為172.009萬 m3/a；邢家渡干渠滲漏系數(shù)為1.010 m2/d·m，側(cè)滲補給總量為52.803萬 m3/a。

(2)根據(jù)區(qū)內(nèi)水位埋深和包氣帶巖性計算本區(qū)潛水蒸發(fā)總量為3 650.671萬 m3/a；區(qū)內(nèi)地下水開采量為775.818萬 m3/a；區(qū)內(nèi)主要排泄河流為小清河，經(jīng)計算，本區(qū)地下水排泄小清河總量為48.681萬 m3/a。地下水側(cè)向徑流排泄總量為137.547萬 m3/a。

2.2 數(shù)學(xué)模型

對于非均質(zhì)各向同性的三維非穩(wěn)定流地下水模型，其數(shù)學(xué)方程描述為：

(1)

定解條件：

H(x，y，z)|t=0=H0(x，y，z)，(x，y，z)∈Ω

(2)

H(x，y，z，t)|(x，y，z)∈B1=H1(x，y，z，t)，(x，y，z)∈B1，t>0

(3)

(4)

式中：Kxx，Kyy，Kzz為xyz方向的滲透系數(shù)，LT-1；h為水頭，L；w為單位體積含水層在單位時間流出或流入地下水的體積，T-1；Ss為儲水系數(shù)；Ω為初始條件，即給定邊界所有點的定水頭；B1為第一類邊界條件；B2為第二類邊界條件。

3 地下水數(shù)值模型

3.1 模型的建立

采用GMS中MODFLOW模塊建立歷城區(qū)地下水流數(shù)值模型，將研究區(qū)空間離散為150×150個單元格，其中有效單元格19 726個。根據(jù)已有資料確定2018年6月-2020年6月為模擬期，時間步長為30 d，共24個應(yīng)力期。導(dǎo)入高程數(shù)據(jù)，輸入各源匯項數(shù)據(jù)。研究區(qū)網(wǎng)格剖分如圖3。

圖3 研究區(qū)三維網(wǎng)格剖分圖

根據(jù)研究區(qū)的地質(zhì)構(gòu)造、巖性及富水性將研究區(qū)按水文地質(zhì)參數(shù)進行分區(qū)[8,9]，如圖4，初始滲透系數(shù)和給水度取值見表1。

圖4 研究區(qū)參數(shù)分區(qū)圖

表1 水文地質(zhì)參數(shù)初值表

根據(jù)研究區(qū)內(nèi)觀測水位進行插值擬合出2018年6月淺層孔隙水流場作為初始流場，如圖5。

圖5 研究區(qū)初始流場擬合圖

3.2 模型識別檢驗

傳統(tǒng)的水文地質(zhì)參數(shù)分區(qū)存在較強的主觀性，忽略了滲透系數(shù)的空間差異性以達到局部最優(yōu)解，而PEST模塊在初始參數(shù)分區(qū)的基礎(chǔ)上可通過添加散點反演滲透系數(shù)場從而達到全局最優(yōu)解。

選取2019年1月至2019年12月觀測井的水位數(shù)據(jù)進行PEST參數(shù)反演，所得滲透系數(shù)場如圖6，觀測井空間分布如圖7，并以此進行地下水動態(tài)檢驗。典型長觀井水位的實測值與計算值擬合情況如圖8，可以看出實測值與模擬計算值動態(tài)變化趨勢一致，擬合情況較好，可以此模型進行地下水動態(tài)預(yù)測[10]。

圖6 滲透系數(shù)反演場 圖7 研究區(qū)觀測井空間分布圖

圖8 研究區(qū)典型長觀井水位過程線擬合圖

4 地下水動態(tài)及流場預(yù)測

地下水開采應(yīng)以可持續(xù)發(fā)展的觀點，在兼顧社會經(jīng)濟和環(huán)境保護的基礎(chǔ)上，合理有效的開發(fā)利用水資源。因此，在已有資料的基礎(chǔ)上，設(shè)計不同的開采方案，通過運行地下水數(shù)值模型來預(yù)測不同開采方案下的地下水流場。

4.1 方案設(shè)計

從充分利用地下水資源的角度考慮，在現(xiàn)狀開采量2.13萬 m3/d的基礎(chǔ)上分別增采10%，20%，50%，70%，預(yù)測5 a后的淺層地下水流場。預(yù)測起始時間為2020年6月，目標時間為2025年6月。預(yù)測期間，大氣降水按多年平均降水量647.9 mm取值。

4.2 結(jié)果分析

(1)經(jīng)過模型預(yù)測，在現(xiàn)狀開采量的基礎(chǔ)上增采10%的情況下，5年后淺層地下水流場見圖9，研究區(qū)內(nèi)典型長觀孔水位動態(tài)變化見圖13(a)。可以看出，當設(shè)計開采量為2.343萬 m3/d時，地下水流場水位變幅在0.1～0.4 m之內(nèi)，地下水降落漏斗區(qū)域面積減小，水位有明顯抬升。

圖9 增采10%地下水流場

(2)在現(xiàn)狀開采量的基礎(chǔ)上增采20%的情況下，5年后淺層地下水流場見圖10，研究區(qū)內(nèi)典型長觀孔水位動態(tài)變化見圖13(b)。相比于增采10%，研究區(qū)內(nèi)典型長觀孔水位抬升變緩，水位變幅在0.1～0.4 m之間。地下水降落漏斗區(qū)域面積仍在減小。

圖10 增采20%地下水流場

(3)在現(xiàn)狀開采量的基礎(chǔ)上增采50%的情況下，5年后淺層地下水流場見圖11，研究區(qū)內(nèi)典型長觀孔水位動態(tài)變化見圖13(c)。當設(shè)計開采量達到3.195萬 m3/d時，研究區(qū)內(nèi)典型長觀孔水位于開采4 a后出現(xiàn)水位下降的趨勢，地下水流場水位變幅處于0.05～0.35 m之間。

圖11 增采50%地下水流場

(4)在現(xiàn)狀開采量的基礎(chǔ)上增采70%的情況下，5年后淺層地下水流場見圖12，研究區(qū)內(nèi)典型長觀孔水位動態(tài)變化見圖13(d)?？梢钥闯?，當設(shè)計增采量達到1.49萬 m3/d時，研究區(qū)內(nèi)典型長觀孔水位于開采2 a后出現(xiàn)水位下降的情況，且下降趨勢較增采50%更為明顯。地下水流場水位變幅仍處于0.05～0.35 m之間。地下水降落漏斗中心區(qū)域開始出現(xiàn)水位下降的情況，可以預(yù)見，若持續(xù)以設(shè)計70%增采量開采淺層地下水將會出現(xiàn)降落漏斗面積擴大，研究區(qū)水位持續(xù)下降等環(huán)境問題。

圖12 增采70%地下水流場

圖13 不同開采方案下典型長觀孔未來5a地下水動態(tài)變化曲線

5 結(jié)語

通過數(shù)值模型預(yù)測，將歷城區(qū)淺層孔隙水開采量控制在3.195萬 m3/d之內(nèi)，可以有效地控制研究區(qū)淺層地下水水位下降，防止地下水降落漏斗持續(xù)擴大，進而引起地面沉降和地下水污染等環(huán)境地質(zhì)問題，達到地下水資源充分且可持續(xù)利用的目標。