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

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

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

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

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

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

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

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

2 地下水流模型

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

2.1.1 地層概化

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

2.1.2 邊界條件

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

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

2.1.3 源匯項(xiàng)

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

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

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

對(duì)于非均質(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為單位體積含水層在單位時(shí)間流出或流入地下水的體積，T-1；Ss為儲(chǔ)水系數(shù)；Ω為初始條件，即給定邊界所有點(diǎn)的定水頭；B1為第一類(lèi)邊界條件；B2為第二類(lèi)邊界條件。

3 地下水?dāng)?shù)值模型

3.1 模型的建立

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

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

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

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

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

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

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

3.2 模型識(shí)別檢驗(yàn)

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

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

圖6 滲透系數(shù)反演場(chǎng) 圖7 研究區(qū)觀(guān)測(cè)井空間分布圖

圖8 研究區(qū)典型長(zhǎng)觀(guān)井水位過(guò)程線(xiàn)擬合圖

4 地下水動(dòng)態(tài)及流場(chǎng)預(yù)測(cè)

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

4.1 方案設(shè)計(jì)

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

4.2 結(jié)果分析

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

圖9 增采10%地下水流場(chǎng)

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

圖10 增采20%地下水流場(chǎng)

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

圖11 增采50%地下水流場(chǎng)

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

圖12 增采70%地下水流場(chǎng)

圖13 不同開(kāi)采方案下典型長(zhǎng)觀(guān)孔未來(lái)5a地下水動(dòng)態(tài)變化曲線(xiàn)

5 結(jié)語(yǔ)

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