濟南市歷城區地下水數值模擬與預測

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

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

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

1 研究區地質概況

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

圖1 研究區地理位置圖

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

圖2 黃河沖積平原區水文地質剖面圖

2 地下水流模型

2.1 水文地質概念模型

2.1.1 地層概化

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

2.1.2 邊界條件

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

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

2.1.3 源匯項

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

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

2.2 數學模型

對于非均質各向同性的三維非穩定流地下水模型，其數學方程描述為：

(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方向的滲透系數，LT-1；h為水頭，L；w為單位體積含水層在單位時間流出或流入地下水的體積，T-1；Ss為儲水系數；Ω為初始條件，即給定邊界所有點的定水頭；B1為第一類邊界條件；B2為第二類邊界條件。

3 地下水數值模型

3.1 模型的建立

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

圖3 研究區三維網格剖分圖

根據研究區的地質構造、巖性及富水性將研究區按水文地質參數進行分區[8,9]，如圖4，初始滲透系數和給水度取值見表1。

圖4 研究區參數分區圖

表1 水文地質參數初值表

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

圖5 研究區初始流場擬合圖

3.2 模型識別檢驗

傳統的水文地質參數分區存在較強的主觀性，忽略了滲透系數的空間差異性以達到局部最優解，而PEST模塊在初始參數分區的基礎上可通過添加散點反演滲透系數場從而達到全局最優解。

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

圖6 滲透系數反演場 圖7 研究區觀測井空間分布圖

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

4 地下水動態及流場預測

地下水開采應以可持續發展的觀點，在兼顧社會經濟和環境保護的基礎上，合理有效的開發利用水資源。因此，在已有資料的基礎上，設計不同的開采方案，通過運行地下水數值模型來預測不同開采方案下的地下水流場。

4.1 方案設計

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

4.2 結果分析

(1)經過模型預測，在現狀開采量的基礎上增采10%的情況下，5年后淺層地下水流場見圖9，研究區內典型長觀孔水位動態變化見圖13(a)?？梢钥闯觯斣O計開采量為2.343萬 m3/d時，地下水流場水位變幅在0.1～0.4 m之內，地下水降落漏斗區域面積減小，水位有明顯抬升。

圖9 增采10%地下水流場

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

圖10 增采20%地下水流場

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

圖11 增采50%地下水流場

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

圖12 增采70%地下水流場

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

5 結語

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