基于MODFLOW與氣候模式的礦區(qū)地下水流模擬

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

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

0 引 言

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

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

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

1 研究區(qū)概況

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

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

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ù)對(duì)地層資料及鉆孔資料的分析，將模型分為七層。進(jìn)行三維地質(zhì)建模時(shí)首先利用CAD 提取原始地形圖的地表高程數(shù)據(jù)，將帶有高程的地形點(diǎn)坐標(biāo)輸入到建模軟件GEO5 中自動(dòng)生成三維地形模型；其次輸入鉆孔資料并對(duì)各個(gè)巖層的顏色及名稱進(jìn)行標(biāo)注，創(chuàng)建水平地層模型；最后生成如圖1所示的三維地質(zhì)模型。綜合地形、地勢(shì)、地質(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 煤層，位于潛水含水層下部的隔水層中，煤炭開采主要會(huì)對(duì)潛水含水層產(chǎn)生影響，因此，本次模擬含水層為第四系松散層孔隙潛水含水層。潛水含水層巖性主要為風(fēng)積砂、灰黃色、淺黃色黃土、殘坡積砂土及沖洪積砂礫石層，厚度為8～63.2 m；下部隔水層為灰白色、褐紫色泥巖、砂質(zhì)泥巖及粉砂巖，厚度為20～50 m。含水層結(jié)構(gòu)示意圖如圖2所示。

2.1.3 邊界條件的概化

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

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

2.1.4 水力特征概化

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

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

對(duì)于非均質(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為含水層垂直補(bǔ)給強(qiáng)度，m/d；E為地下水蒸發(fā)排泄強(qiáng)度，m/d；P為含水層開采強(qiáng)度，m/d；h0為含水層的初始水頭，m；q為含水層側(cè)向單寬補(bǔ)排量，m2/d，流入為正，流出為負(fù)，隔水邊界為0；Γ1為隔水邊界；Γ2為流量邊界；n→為邊界上的外法線方向［19］。

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

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

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

2.4 模型識(shí)別及驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 煤炭開采對(duì)地下水水位的影響分析

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

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

模擬結(jié)果顯示，隨著多年對(duì)研究區(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)狀地下水補(bǔ)給總量為573.22 萬m3，排泄總量為1 038.5 萬m3，均衡差為465.28 萬m3；2035年地下水補(bǔ)給總量為354.18 萬m3，排泄總量為419.09 萬m3，均衡差為64.91 萬m3。通過表4可以看出，同2018年相比，2035年的各項(xiàng)水量均減少，且補(bǔ)給量與排泄量的均衡差也明顯減小，表明研究區(qū)的地下水系統(tǒng)逐漸趨向于新的平衡。一般情況下隨著地下水位的降低會(huì)激發(fā)周邊地下水的側(cè)向補(bǔ)給，而本文的側(cè)向流入?yún)s大幅度減少，這是由于多年的過度開采地下水，導(dǎo)致研究區(qū)周邊地下水也相應(yīng)地受到影響，進(jìn)而無法滿足研究區(qū)地下水系統(tǒng)的側(cè)向補(bǔ)給。

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

4 結(jié) 論

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

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

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