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

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

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

0 引 言

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

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

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

1 研究區概況

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

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

2 地下水流模型構建

2.1 水文地質概念模型

2.1.1 三維地質模型構建

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

2.1.2 含水層結構概化

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

2.1.3 邊界條件的概化

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

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

2.1.4 水力特征概化

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

2.2 數學模型

對于非均質、各向異性、空間三維結構、非穩定地下水流系統，數學模型如下［18］：

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

2.3 模型的離散及參數分區

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

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

2.4 模型識別及驗證

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

表2 研究區水文地質參數表Tab.2 Hydrogeological parameters of the study area

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

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

3 未來氣候變化情景構建

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

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

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

3.1 模型預測

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

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

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

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

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

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

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

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

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

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

4 結 論

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

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

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