許昌市建安區(qū)北部采空區(qū)三維地下水流數(shù)值模擬

王靈敏, 王軍強(qiáng), 劉榮慧, 李卓然, 常偉麗， 鄭萍

(河南省地質(zhì)礦產(chǎn)勘查開發(fā)局第二地質(zhì)勘查院,河南 鄭州 450045)

近年來，黃昕霞等[1]、MARC-ANDRE Lavigne等[2]、ELY D M等[3]、吳樂等[4]、王建紅等[5]先后應(yīng)用數(shù)值模擬技術(shù)分別建立了山西運(yùn)城盆地、紐約沙托河流域、哥倫比亞高原地區(qū)、北京市西山地區(qū)、黑龍江黑河流域中游平原區(qū)的地下水流數(shù)值模型，對(duì)其地下水水位變化特征進(jìn)行了模擬研究。筆者利用數(shù)字模型對(duì)地下水位變化特征進(jìn)行研究。

文中研究區(qū)位于河南省許昌市建安區(qū)北部，為非工業(yè)區(qū)，主要人類活動(dòng)為道路建設(shè)、農(nóng)業(yè)生產(chǎn)、礦山開采等。蘇橋鎮(zhèn)武莊村南部許昌鐵礦于2013年初開始開采，在許昌鐵礦礦山基本設(shè)施建設(shè)期間,建安區(qū)北部出現(xiàn)地下水水位大幅下降、地面塌陷、地面沉降等環(huán)境地質(zhì)災(zāi)害問題。2013—2015年，該研究區(qū)地下水水位降幅達(dá)2～3 m；2013年,武莊村發(fā)生地面塌陷，并造成1人死亡；2015年,礦區(qū)內(nèi)水泥道路突然發(fā)生崩裂錯(cuò)斷，形成9處陷洞，部分房屋墻體產(chǎn)生裂縫[6]。上述環(huán)境地質(zhì)問題影響了當(dāng)?shù)鼐用竦纳a(chǎn)和生活。

該地區(qū)盛產(chǎn)鐵礦，目前針對(duì)這些鐵礦開展的大多數(shù)工作集中在礦產(chǎn)勘查及環(huán)境地質(zhì)問題初步調(diào)查階段，對(duì)當(dāng)?shù)丨h(huán)境地質(zhì)問題的發(fā)展趨勢(shì)及原因尚未有研究。本文在前人調(diào)查資料的基礎(chǔ)上,應(yīng)用當(dāng)代水文地質(zhì)學(xué)研究和應(yīng)用的最重要手段之一——地下水?dāng)?shù)值模擬技術(shù)[7]對(duì)研究區(qū)內(nèi)地下水流問題進(jìn)行了研究，建立了研究區(qū)的三維水文地質(zhì)模型，并對(duì)該區(qū)域地下水水位動(dòng)態(tài)及變化趨勢(shì)進(jìn)行了研究。研究結(jié)果將為該研究區(qū)合理開采地下水資源提供參考。

1 研究區(qū)的水文地質(zhì)條件

研究區(qū)為山前沖洪積平原，地勢(shì)平緩，西北向東南地形坡度為5‰。研究區(qū)地下水的賦存和分布規(guī)律以地層巖性為基礎(chǔ)[8]，區(qū)內(nèi)廣泛沉積新生界巨厚松散沉積層，厚度50～400 m。松散孔隙含水層劃分為3個(gè)含水段(In1、In2、In3)，如圖1所示。

圖1 研究區(qū)水文地質(zhì)剖面圖

圖1中自上而下劃分為：第四系孔隙潛水、淺層承壓水含水段(In1)、新近系孔隙承壓含水段(In2)和古近系孔隙承壓含水段(In3)[9]。含水段(In1)，含水巖組為第四系鈣質(zhì)結(jié)核層，含鈣質(zhì)結(jié)核粉土，為中等富水區(qū)，單井5 m降深涌水量100～500 m3/d。含水段(In2、In3)，含水巖組為沖洪積粉細(xì)砂、中砂層。研究區(qū)內(nèi)含水段In2和In3為富水區(qū)，單井15 m降深涌水量1 000～3 000 m3/d[10]。各含水段間隔水層巖性主要為粉質(zhì)黏土、弱成巖粉質(zhì)黏土，隔水性較好。古近系孔隙承壓含水段In3底部為古近系泥礫巖和泥質(zhì)角礫巖，與基巖呈不整合接觸，隔水性較好。本次模擬對(duì)象為松散孔隙含水層。新近系含水段In2與古近系含水段In3均為承壓含水層，地下水補(bǔ)給來源與排泄方式一致，因此，模擬中將新近系含水段In2與古近系含水段In3合并為一個(gè)含水段進(jìn)行研究。

研究區(qū)內(nèi)第四系孔隙潛水、淺層承壓水含水段(In1)主要接受大氣降水補(bǔ)給和少量徑流補(bǔ)給；以蒸發(fā)、農(nóng)田供水及礦區(qū)排水形式排泄。新近系孔隙承壓含水段(In2)和古近系孔隙承壓含水段(In3)主要接受地下水側(cè)滲補(bǔ)給，以礦區(qū)排水及地下徑流的方式排泄。研究區(qū)內(nèi)地下水流方向?yàn)槲鞅薄獤|南向。

2 研究區(qū)地下水的三維數(shù)學(xué)模型

在分析研究區(qū)水文地質(zhì)條件的基礎(chǔ)上，應(yīng)用連續(xù)性原理及達(dá)西定律，構(gòu)建了研究區(qū)地下水的三維非穩(wěn)定流數(shù)學(xué)模型[11-13]，即：

(1)

3 研究區(qū)地下水的數(shù)值模擬

使用地下水模型軟件GMS(Groundwater Modeling System)中的Modflow模塊建立地下水流數(shù)值模型。GMS是目前國際上通用的地下水模擬軟件，可用來模擬飽和、非飽和流條件下的地下水流和溶質(zhì)運(yùn)移等情況[14]。Modflow作為地下水流數(shù)值模擬軟件，實(shí)現(xiàn)了概念模型和數(shù)值模型的無縫轉(zhuǎn)化，具備強(qiáng)大的二維和三維可視化功能，可以很好地模擬水井、河流、溪流等對(duì)非均質(zhì)地質(zhì)條件和復(fù)雜邊界條件的水位系統(tǒng)的影響[15]。

模擬區(qū)邊界：東以武莊礦區(qū)東部邊界為界，北以石固村南部邊界為界，西以武莊礦區(qū)西部邊界為界，南以石梁河為界，總面積為9.1 km2。邊界條件概化如圖2所示：北部、西部為模型邊界1，定為流入邊界；東部為模型邊界2，定為流出邊界；南部為模型邊界3，定為定水頭邊界。

根據(jù)研究區(qū)的含水層結(jié)構(gòu)、邊界條件和地下水流場(chǎng)特征，垂向上將地下水系統(tǒng)結(jié)構(gòu)概化為3層，如圖3所示。圖3中：第一層為第四系松散潛水—微承壓含水層，第二層為黏土隔水層，第三層為新近系松散土層承壓含水層。模擬區(qū)剖分為39行、51列規(guī)則網(wǎng)格，3層，各層均采用118 m×118 m的剖分格式。其中有效單元為4 437個(gè)，如圖4所示。模型的模擬期為2015年9月—2021年9月，識(shí)別驗(yàn)證期為2015年9月—2018年9月，模擬預(yù)測(cè)期為2018年9月—2021年9月。將整個(gè)預(yù)測(cè)期劃分為3個(gè)應(yīng)力期，每個(gè)應(yīng)力期為一個(gè)相應(yīng)的自然月，計(jì)算的時(shí)間步長(zhǎng)為15 d。模擬初始水位采用2015年9月的實(shí)測(cè)水位(只做一次水位測(cè)量)，識(shí)別驗(yàn)證采用2018年9月的(只做一次水位測(cè)量)實(shí)測(cè)水位。

圖3 三維地下水概念模型界面圖

圖4 研究區(qū)單元剖分網(wǎng)格圖

模型采用“試錯(cuò)法”(trial-and-error)來間接率定水流模型的相關(guān)參數(shù)[16]。對(duì)含水層和弱透水層進(jìn)行參數(shù)分區(qū)，給定初始值，具體見表1。圖5為研究區(qū)2015年9月的實(shí)測(cè)水位,即研究區(qū)地下水?dāng)?shù)值模擬的初始水位。在識(shí)別驗(yàn)證期，根據(jù)實(shí)際水位與模擬計(jì)算水位的比對(duì)結(jié)果，不斷調(diào)整參數(shù)，直到模擬水位與實(shí)際水位的誤差達(dá)到允許范圍以內(nèi)。研究區(qū)有13眼觀測(cè)井進(jìn)行水位擬合，觀測(cè)井分布位置如圖6所示，基本覆蓋研究區(qū)。

表1 含水層及隔水層參數(shù)取值范圍

圖5 研究區(qū)初始水位(單位：m)

圖6 研究區(qū)地下水觀測(cè)井分布圖

在識(shí)別驗(yàn)證期末(即2018年9月)，對(duì)研究區(qū)13眼觀測(cè)井的模型計(jì)算水位值與實(shí)際觀測(cè)水位值進(jìn)行對(duì)比擬合。識(shí)別驗(yàn)證期末研究區(qū)觀測(cè)井實(shí)際觀測(cè)水位值與模型計(jì)算水位值見表2。由表2知，二者水位值擬合良好，地下水水位標(biāo)高誤差小于0.4 m。

除此之外，在識(shí)別驗(yàn)證期末，對(duì)研究區(qū)模擬流場(chǎng)與實(shí)際流場(chǎng)進(jìn)行了擬合，結(jié)果如圖7所示。由圖7可知，模擬流場(chǎng)與實(shí)際流場(chǎng)總體流動(dòng)方向相同，大部分節(jié)點(diǎn)水位擬合較好，模型參數(shù)及源匯項(xiàng)設(shè)置符合實(shí)際情況。

識(shí)別驗(yàn)證期內(nèi)研究區(qū)的地下水水位變幅如圖8所示。其中，2018年9月的地下水水位為模型計(jì)算水位，2015年9月的水位為實(shí)際觀測(cè)水位。由圖8可知，識(shí)別驗(yàn)證期內(nèi)，研究區(qū)地下水水位整體變幅較小，研究區(qū)東部部分地區(qū)如磨李、西張、雷莊區(qū)域的地下水水位呈0.0～0.5 m的小幅下降，武莊塌陷區(qū)的地下水水位呈0.0～0.5 m的小幅上升，其他區(qū)域的地下水水位呈0.5～1.0 m的幅度上升。

表2 研究區(qū)觀測(cè)井2018年的觀測(cè)水位與模型計(jì)算水位 m

圖7 研究區(qū)水位擬合圖(單位：m)

圖8 2018年9月的預(yù)測(cè)水位較2015年9月的實(shí)測(cè)水位的變幅圖

識(shí)別驗(yàn)證期內(nèi)研究區(qū)地下水水位變幅較小的原因?yàn)椋孩俑鶕?jù)當(dāng)?shù)貧庀缶仲Y料，研究區(qū)2015—2018年的降雨量較2013—2015年的降雨量有所增加，增加量為257 mm/a，因此，研究區(qū)地下水的降雨入滲補(bǔ)給量較2013—2015年的降雨入滲補(bǔ)給量大幅增加；②研究區(qū)南部石梁河由于礦井排水，水位較高，常年滲漏補(bǔ)給地下水；③由于模擬區(qū)屬于山前徑流區(qū)，再加上7～9月雨季入滲的補(bǔ)給，西部水位上升幅度大，水力梯度增大，徑流排泄增強(qiáng)，因此，東部地區(qū)水位有小幅下降。

4 模型預(yù)測(cè)

建立水流模型的目的是為了預(yù)測(cè)在設(shè)定地下水開采條件下研究區(qū)內(nèi)地下水流的動(dòng)態(tài)變化特征，從而判斷該種開采方案是否科學(xué)，為合理利用地下水提供建議。2012—2015年,研究區(qū)出現(xiàn)水位大幅下降、地面塌陷及地裂縫等一系列環(huán)境地質(zhì)問題。目前，研究區(qū)內(nèi)水位下降的原因及地下水流動(dòng)態(tài)特征亟待解決。模型預(yù)測(cè)期選為2018—2021年，開采方案為現(xiàn)狀開采，即保持現(xiàn)有2018年的礦井抽水量和地下水開采項(xiàng)、開采量、開采井位置不變。

2021年9月與2018年9月地下水水位對(duì)比如圖9所示。由圖9可知，在礦井保持2018年地下水抽水量不變的條件下，武莊區(qū)域地下水水位降落漏斗區(qū)的中心水位出現(xiàn)上升，該地下水水位降落漏斗區(qū)域的地下水水位升幅最大，為1.0～1.5 m；雷莊區(qū)域地下水水位降落漏斗區(qū)的中心水位出現(xiàn)上升，該地下水水位降落漏斗區(qū)域的地下水水位升幅為0.5～1.0 m；研究區(qū)其它大部分區(qū)域的地下水水位也趨于平緩上升，上升幅度小于0.5 m,地下水水位的變化幅度較小。說明研究區(qū)地下水水位相對(duì)穩(wěn)定。

圖9 2021年9月的預(yù)測(cè)水位較2018年9月的實(shí)測(cè)水位的變幅圖

上述表明，2013—2015年，武莊區(qū)的礦井運(yùn)營初期，因礦井大量抽水，導(dǎo)致研究區(qū)內(nèi)地下水排泄量大于補(bǔ)給量，研究區(qū)內(nèi)地下水水位大幅下降；2015—2018年，研究區(qū)地下水形成地下水水位降落漏斗，激發(fā)研究區(qū)外圍補(bǔ)水，補(bǔ)給量大于排泄量，地下水水位基本穩(wěn)定；預(yù)測(cè)期的2018—2021年，地下水補(bǔ)給量依然大于地下水排泄量，地下水水位降落漏斗區(qū)的中心水位上升，區(qū)域地下水水位呈整體小幅上升趨勢(shì)。

5 結(jié)語

1)通過資料收集、水位統(tǒng)調(diào)、抽水試驗(yàn)等水文地質(zhì)調(diào)查工作，應(yīng)用GMS(Modflow模塊)地下水流模擬軟件建立了研究區(qū)松散孔隙含水層三維數(shù)值模型，模型經(jīng)過識(shí)別驗(yàn)證，能夠較好地模擬研究區(qū)的地下水流。

2)模型模擬結(jié)果顯示，2015—2018年，研究區(qū)的地下水水位變幅較小，武莊礦區(qū)周圍區(qū)域的地下水水位升幅在0.5 m左右，模擬區(qū)東部雷莊區(qū)域的地下水水位降落漏斗區(qū)與西張區(qū)域的地下水水位降落漏斗區(qū)，這兩個(gè)地區(qū)的地下水水位有0.2～0.5 m的下降。

3)模擬預(yù)測(cè)結(jié)果顯示，在現(xiàn)狀開采條件下，即礦井保持2018年抽水量繼續(xù)抽水，研究區(qū)大部分地區(qū)的地下水水位趨于穩(wěn)定及平緩上升，上升幅度小于0.5 m。武莊區(qū)域的地下水水位降落漏斗區(qū)的中心水位升高，區(qū)域地下水水位升幅為1.0～1.5 m；研究區(qū)東北部雷莊區(qū)域的地下水水位降落漏斗區(qū)的水位也略有回升，升幅為0.5～1.0 m。

4)2012—2013年，研究區(qū)由于礦井抽水形成的地下水水位降落漏斗已經(jīng)激發(fā)了遠(yuǎn)程的側(cè)向補(bǔ)給，研究區(qū)內(nèi)地下水水位降落漏斗中心點(diǎn)水位升高，區(qū)域水位整體上升，水位降落漏斗進(jìn)入恢復(fù)涵養(yǎng)期[17]，研究區(qū)地下水已形成新的均衡。