劉中培，李 鑫，陳 瑩，韓宇平

（1.華北水利水電大學(xué)，鄭州450045；2.煤炭開采水資源保護與利用國家重點實驗室，北京100011；3.水利部水資源管理中心，北京100038）

0 引 言

煤炭是重要的礦產(chǎn)資源，由于我國國民經(jīng)濟近年來取得高速發(fā)展，煤炭的需求量也隨之與日俱增［1］。為了保障煤礦的正常開采以及礦區(qū)的生活用水需求，要對煤礦進行疏干排水，以解決煤礦開采所引起的涌水問題［2］。長時間的疏干排水會導(dǎo)致礦區(qū)周圍地下水水位下降，打破地下水的平衡狀態(tài)，破壞地下水資源，進一步加劇地下水短缺的局面［3］。因此，研究煤礦開采對地下水水位的影響對于地下水資源的可持續(xù)利用具有重大意義［4］。

國內(nèi)外的眾多學(xué)者就煤炭開采對地下水水位影響的問題做了大量研究。陳社明等［5］通過Visual MODFLOW 軟件構(gòu)建了研究區(qū)的地下水流數(shù)值模型，并以此預(yù)測了未來研究區(qū)各含水層水位在煤炭開采過程中的變化特征，分析礦山排水對不同含水層水位的影響。BISWAS H 等［6］在分析沃克爾盆地水文地質(zhì)條件和水資源利用的基礎(chǔ)上利用MODFLOW 構(gòu)建了研究區(qū)地下水流數(shù)值模型，并利用模型模擬了3 種場景來研究補給減少和開采量增加對地下水水位的影響。白曉等［7］將ARIMA 模型預(yù)測的降水量用于MODFLOW 的預(yù)測模型之中，對研究區(qū)巖溶地下水資源量和水位動態(tài)變化進行了模擬和預(yù)測。汪麗芳等［8］建立了多種開采情境下的地下水位動態(tài)預(yù)測模型，并利用多元回歸方法分析地下水水位與地下水開采量之間的關(guān)系。這些研究成果為地下水?dāng)?shù)值模擬技術(shù)的應(yīng)用提供了重要的參考依據(jù)，但是傳統(tǒng)的利用數(shù)值模擬技術(shù)預(yù)測地下水水位往往忽略了未來氣候變化及煤礦的采掘進度，預(yù)測結(jié)果不能較準(zhǔn)確地反映地下水水位隨著氣候及采掘進度的動態(tài)變化［9］。

本文根據(jù)研究區(qū)的地下水水位監(jiān)測數(shù)據(jù)和詳細的水文地質(zhì)資料，利用Visual MODFLOW 建立研究區(qū)地下水流數(shù)值模型，并依靠實測數(shù)據(jù)對模型進行識別和驗證［10］。最后利用驗證好的模型結(jié)合氣候模式所預(yù)測的未來降水量以及研究區(qū)的采煤規(guī)劃，預(yù)測研究區(qū)在煤炭開采條件下潛水含水層的水位變化趨勢，探究煤炭開采對地下水水位的影響。

1 研究區(qū)概況

研究區(qū)位于鄂爾多斯市伊金霍洛旗烏蘭木倫鎮(zhèn)境內(nèi)，包括烏蘭木倫煤礦和柳塔煤礦，屬于干旱半干旱大陸性氣候，四季寒暑巨變。研究區(qū)位于毛烏素沙漠的邊緣地帶，陜北黃土高原的東北部，地勢北高南低，東高西低。從總體趨勢看，研究區(qū)發(fā)育的小構(gòu)造受褶皺構(gòu)造控制明顯，一般發(fā)育于褶皺的翼部，基本呈北西與北東兩組。區(qū)內(nèi)地表水系發(fā)育，西北緊鄰公涅爾蓋溝，西南倚靠烏蘭木倫河。

研究區(qū)內(nèi)地下水可分為新生界松散層孔隙潛水和中生界碎屑巖類孔隙、裂隙水兩大類。松散層孔隙潛水含水層是區(qū)內(nèi)主要的含水層，包括地勢較高處的風(fēng)積砂含水層和溝谷地帶的沖洪積湖積砂礫石層含水層，接受大氣降水和地表水的補給，沿基巖頂面流動，在溝谷處排泄。基巖裂隙潛水——承壓水含水層主要為侏羅系地層中的砂巖，潛水區(qū)補給來源為大氣降水的直接補給和松散層孔隙潛水含水層的補給，沿地層傾向流動；承壓水區(qū)補給來源較為單一，主要為徑流補給和基巖裸露區(qū)段大氣降水的直接補給，徑流排泄、泉水排泄是該含水層的持久排泄方式。

2 地下水流模型構(gòu)建

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

2.1.1 三維地質(zhì)模型構(gòu)建

研究區(qū)地表大部分被第四系覆蓋，僅在烏蘭木倫煤礦西鄰的烏蘭木倫河?xùn)|岸有少量基巖出露。區(qū)內(nèi)沉積地層由老到新為：三疊系上統(tǒng)延長組、侏羅系中下統(tǒng)延安組、侏羅系中統(tǒng)直羅組、侏羅系中統(tǒng)安定組、白堊系下統(tǒng)伊金霍洛組、第四系上更新統(tǒng)、第四系全新統(tǒng)。根據(jù)對地層資料及鉆孔資料的分析，將模型分為七層。進行三維地質(zhì)建模時首先利用CAD 提取原始地形圖的地表高程數(shù)據(jù)，將帶有高程的地形點坐標(biāo)輸入到建模軟件GEO5 中自動生成三維地形模型；其次輸入鉆孔資料并對各個巖層的顏色及名稱進行標(biāo)注，創(chuàng)建水平地層模型；最后生成如圖1所示的三維地質(zhì)模型。綜合地形、地勢、地質(zhì)條件等構(gòu)建的三維地質(zhì)模型基本展現(xiàn)了礦區(qū)的地質(zhì)特征［11］。并且此模型通過任意切換視角和剖切等手段可以清晰、真切地反映地層復(fù)雜的內(nèi)部結(jié)構(gòu)，為水文地質(zhì)模型的構(gòu)建奠定了基礎(chǔ)［12］。

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

研究區(qū)主采煤層為1～2 煤層，位于潛水含水層下部的隔水層中，煤炭開采主要會對潛水含水層產(chǎn)生影響，因此，本次模擬含水層為第四系松散層孔隙潛水含水層。潛水含水層巖性主要為風(fēng)積砂、灰黃色、淺黃色黃土、殘坡積砂土及沖洪積砂礫石層，厚度為8～63.2 m；下部隔水層為灰白色、褐紫色泥巖、砂質(zhì)泥巖及粉砂巖，厚度為20～50 m。含水層結(jié)構(gòu)示意圖如圖2所示。

2.1.3 邊界條件的概化

垂向邊界的概化：因為潛水含水層與外界發(fā)生降雨入滲、蒸發(fā)排泄等垂向水量交換，所以模擬區(qū)頂部以潛水含水層的自由水面為系統(tǒng)上邊界，由于研究區(qū)的最低潛水水位埋深仍為12.7 m，故本次模擬中忽略了蒸散發(fā)等帶來的影響；底部以泥巖、砂質(zhì)泥巖等構(gòu)成的隔水層為系統(tǒng)下邊界［13］。

側(cè)向邊界的概化：研究區(qū)西北部的公涅爾蓋溝上游常年有水，對潛水含水層存在補給，概化為流量邊界（流入邊界），公涅爾蓋溝下游地勢較低，水位低于區(qū)域地下水位，概化為流量邊界（流出邊界）。研究區(qū)西南部的烏蘭木倫河屬于季節(jié)性河流，河流水位常年低于地下水位，概化為流量邊界（流出邊界）；東側(cè)烏蘭木倫煤礦井田邊界保留煤柱，并且位于地表分水嶺附近，概化為隔水邊界。水文地質(zhì)概念模型如圖3所示。

2.1.4 水力特征概化

研究區(qū)的含水層的分布范圍較廣，在常溫常壓下的地下水系統(tǒng)遵循質(zhì)量守恒定律和達西定律［14］。參數(shù)隨著空間的變化而發(fā)生改變，垂直與水平方向上的滲透系數(shù)有所不同，因此將參數(shù)概化為非均質(zhì)各向異性［15］。由于地下水系統(tǒng)的滲流運動要素隨時間和空間變化，故將地下水系統(tǒng)概化為非穩(wěn)定流［16］?？傊?，研究區(qū)可概化為非均質(zhì)、各向異性、三維非穩(wěn)定地下水流系統(tǒng)［17］。

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

對于非均質(zhì)、各向異性、空間三維結(jié)構(gòu)、非穩(wěn)定地下水流系統(tǒng)，數(shù)學(xué)模型如下［18］：

式中：Ω為滲流區(qū)域；H為含水層水頭，m；Z為含水層底板高程，m；K為滲透系數(shù)；Kxx、Kyy、Kzz分別為x、y、z方向上的含水層滲透系數(shù)，m/d；μ*為貯水系數(shù)；W為含水層垂直補給強度，m/d；E為地下水蒸發(fā)排泄強度，m/d；P為含水層開采強度，m/d；h0為含水層的初始水頭，m；q為含水層側(cè)向單寬補排量，m2/d，流入為正，流出為負，隔水邊界為0；Γ1為隔水邊界；Γ2為流量邊界；n→為邊界上的外法線方向［19］。

2.3 模型的離散及參數(shù)分區(qū)

研究區(qū)的面積大小為57.7 km2，在水平方向上將模擬區(qū)域剖分成100 m×100 m的單元格。滲透系數(shù)及給水度等水文地質(zhì)參數(shù)主要通過研究區(qū)水文地質(zhì)試驗和經(jīng)驗值獲得，以鉆孔ZK3509 抽水試驗數(shù)據(jù)為例，詳見表1。根據(jù)試驗數(shù)據(jù)和水文地質(zhì)條件等將含水層共分為7個參數(shù)分區(qū)［20］，如圖4所示。

表1 研究區(qū)典型鉆孔抽水試驗成果表Tab.1 Pumping test results of typical boreholes in the study area

2.4 模型識別及驗證

模型的識別期為2017年1月1日到2017年12月31日。利用2017年的地下水水位實測數(shù)據(jù)進行參數(shù)識別，從水文地質(zhì)參數(shù)的初值出發(fā)，根據(jù)參數(shù)變化范圍和實際水位差值對模型進行分區(qū)反演試算，直到模型的計算值與研究區(qū)的實測值擬合較好為止［3］。最終得到各監(jiān)測井所在點的計算值與實測值的擬合曲線。部分監(jiān)測井的擬合曲線如圖5所示，含水層水文地質(zhì)參數(shù)識別結(jié)果見表2。

表2 研究區(qū)水文地質(zhì)參數(shù)表Tab.2 Hydrogeological parameters of the study area

模型的驗證期為2018年1月1日至2018年5月31日。將2018年實測降水量、開采量等輸入到已經(jīng)建好的模型中，對區(qū)內(nèi)地下水位的計算值和實測值進行比較，如圖6所示，驗證期內(nèi)計算水位接近實測水位，二者相差較小。

模型識別期和驗證期的各監(jiān)測井的計算水位與實測水位的均方誤差分別為0.069 m2和0.015 m2。結(jié)果表明所建模型可以較好地模擬地下水流情況，可應(yīng)用于地下水水位預(yù)測。

3 未來氣候變化情景構(gòu)建

未來氣候預(yù)估是指氣候系統(tǒng)響應(yīng)溫室氣體和氣溶膠的排放、濃度或輻射強迫等情景所作出的預(yù)估，通常是基于氣候模式的模擬結(jié)果［21］。氣候模式是當(dāng)前預(yù)測未來氣候變化情景的重要手段，其模擬數(shù)據(jù)采用ISI-MIP 提供的GFDL-ESM2M，HadGEM2-ES，IPSL-CM5A-LR，MIROC-ESM-CHEM，Nor-ESM1-M 等5 個全球氣候模式日值數(shù)據(jù)的插值和修正的結(jié)果，空間水平分辨率為0.5°×0.5°，時間范圍是1951-2050年［22］。本文選用IPCC AR5 中確定的“典型濃度路徑”（RCPs）作為未來情景，涵蓋了RCP2.6、RCP4.5、RCP6.0和RCP8.5等4個情景［23］。

首先使用實測降水格點數(shù)據(jù)對氣候模式模擬數(shù)據(jù)集模擬的降水進行初步檢驗。實測數(shù)據(jù)來源于1961年1月以來的中國地面水平分辨率0.5° × 0.5°的日值格點數(shù)據(jù)，而氣候模式的基準(zhǔn)期為1951-2000年，所以本文以1961-2000年為基準(zhǔn)期。然后利用中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)提供的實測降水值與氣候模式模擬數(shù)據(jù)集模擬的降水值進行對比，最后通過逐月同比例縮放法對降水值做進一步修正。對于降水結(jié)果的修正，可用如下公式進行求解，見公式（1）和公式（2）。

式中：λGCMs為修正系數(shù)為基準(zhǔn)期內(nèi)逐月多年實測降水平均值為基準(zhǔn)期內(nèi)逐月多年模擬降水平均值；x′GCMs-f，i為未來時段修正后的模擬值；xGCMs-f，i為未來時段逐月模擬值［24］。

3.1 模型預(yù)測

由于地下水水位受氣象因素中降水量的影響較大，此處選用RCP4.5即中等濃度路徑情景作為氣候模式的模擬情景，而氣候模式的選取則要根據(jù)5套氣候模式所模擬的逐月多年平均降水量的相對誤差來確定［14］，計算結(jié)果見表3。

表3 不同氣候模式下研究區(qū)逐月多年平均降水量相對誤差 %Tab.3 Relative error of annual average precipitation in the study area under different climate models

通過表3可以看出由氣候模式GFDL-ESM2M 所計算出的模擬值與實測值最接近，因此研究區(qū)的未來降水量選取RCP4.5情景下GFDL-ESM2M模式中的模擬結(jié)果。

依據(jù)2015-2018年研究區(qū)涌水量實測值以及未來礦區(qū)的采煤規(guī)劃，預(yù)測2020-2035年的地下水開采量。模型預(yù)測的初始水位取模型驗證期的末刻時間，即2018年5月31日的實測水位，利用已經(jīng)建立并驗證的模型預(yù)測未來研究區(qū)的地下水流場變化情況。模型潛水含水層的現(xiàn)狀與未來地下水流場對比，如圖7所示。

通過圖7可以看出在煤礦正常開采的條件下，地下水流場變化較為明顯。研究區(qū)大部分區(qū)域都表現(xiàn)出水位降低的趨勢，疏干面積也隨之增多。東北部區(qū)域出現(xiàn)局部的水位回升，而且水位最高點也向西偏移；西南方向因礦區(qū)開采形成的地下水降落漏斗面積減小，但是降落漏斗的周圍區(qū)域出現(xiàn)大幅度的水位下降。

3.2 煤炭開采對地下水水位的影響分析

礦區(qū)因不斷地開采煤炭而產(chǎn)生的大量礦井涌水不僅影響煤炭的正常開采，還會威脅到人身安全。因此需要及時將涌水排出，但是大量的地下水被集中排出將會導(dǎo)致周邊地下水水位打破其原有的平衡狀態(tài)。為了更加直觀地分析煤炭開采過程中地下水水位的變化情況，對比2018年與2035年的地下水流場圖（圖7）與地下水水量均衡表（表4），分析煤炭開采對潛水水位的影響。

表4 2018與2035年地下水水量均衡 萬m3Tab.4 Groundwater balance in 2019 and 2035

模擬結(jié)果顯示，隨著多年對研究區(qū)地下水集中開采，地下水等水位線由疏變密，地下水流速加快，地下水均衡也發(fā)生變化。從2018-2035年，烏蘭木倫煤礦北部區(qū)域水位回升0.97～17.02 m，東部區(qū)域水位下降11.32～28.27 m，南部區(qū)域下降0.24～28.74 m，西部區(qū)域下降11.32～16.95 m；柳塔煤礦西南區(qū)域的地下水降落漏斗范圍有所減小，漏斗中心水位回升4.03 m，東部及南部區(qū)域水位下降0.24～17.69 m，西部區(qū)域平均下降6.44 m。研究區(qū)現(xiàn)狀地下水補給總量為573.22 萬m3，排泄總量為1 038.5 萬m3，均衡差為465.28 萬m3；2035年地下水補給總量為354.18 萬m3，排泄總量為419.09 萬m3，均衡差為64.91 萬m3。通過表4可以看出，同2018年相比，2035年的各項水量均減少，且補給量與排泄量的均衡差也明顯減小，表明研究區(qū)的地下水系統(tǒng)逐漸趨向于新的平衡。一般情況下隨著地下水位的降低會激發(fā)周邊地下水的側(cè)向補給，而本文的側(cè)向流入?yún)s大幅度減少，這是由于多年的過度開采地下水，導(dǎo)致研究區(qū)周邊地下水也相應(yīng)地受到影響，進而無法滿足研究區(qū)地下水系統(tǒng)的側(cè)向補給。

綜上分析，同現(xiàn)狀情況相比，煤礦經(jīng)多年開采，由于大量疏干排水導(dǎo)致研究區(qū)內(nèi)地下水資源總量嚴重削減，進而造成研究區(qū)內(nèi)的潛水水位普遍下降。隨著礦區(qū)不斷地開采，采空區(qū)增多，并且依據(jù)采煤規(guī)劃，由于采煤工作面的轉(zhuǎn)移和采煤數(shù)量的減少，地下水開采量也隨之減少，故研究區(qū)的局部區(qū)域的潛水水位有所回升。

4 結(jié) 論

（1）根據(jù)翔實的水文地質(zhì)資料及地下水水位監(jiān)測資料所建立的地下水流數(shù)值模型對研究區(qū)地下水水位的計算值與實測值擬合較好，水文地質(zhì)參數(shù)的選取也較為合理，模型可用于地下水流的計算與預(yù)測。將氣候模式所預(yù)測的降水量應(yīng)用于MODFLOW 模型預(yù)測之中，提高了對未來地下水水位變化情況預(yù)測的準(zhǔn)確度。

（2）煤礦在現(xiàn)狀條件下開采，研究區(qū)內(nèi)地下水補給總量為573.22 萬m3，排泄總量為1 038.5 萬m3，均衡差為465.28 萬m3，地下水系統(tǒng)處于負均衡狀態(tài)。由于礦區(qū)長時間的地下水超采，研究區(qū)潛水水位下降嚴重，柳塔煤礦所形成的地下水降落漏斗范圍較大。

（3）模型根據(jù)研究區(qū)的降雨預(yù)測結(jié)果以及未來的采煤規(guī)劃，預(yù)測2020-2035年研究區(qū)的地下水資源總量逐漸減少，潛水水位普遍下降；隨著多年的地下水超采，研究區(qū)多個區(qū)域被疏干，并且隨著采煤工作面的轉(zhuǎn)移以及地下水開采量逐漸減小，局部區(qū)域出現(xiàn)水位回升現(xiàn)象，地下水降落漏斗范圍有所減小?！?/p>