左偉，劉運濤，李建林，王樹威

（1.河南省地質(zhì)礦產(chǎn)勘查開發(fā)局 第五地質(zhì)勘查院，河南 鄭州 450001；2.河南理工大學(xué) 資源環(huán)境學(xué)院，河南 焦作 454000）

0 引言

地?zé)豳Y源是一種清潔、低碳的可再生能源。地?zé)豳Y源的開發(fā)利用可有效優(yōu)化能源供給結(jié)構(gòu)，對解決全球化的生態(tài)環(huán)境問題具有深遠(yuǎn)影響。按照埋藏深度，可將地?zé)豳Y源劃分為淺層地?zé)?、水熱型地?zé)岷透蔁釒r［1］。淺層地?zé)嵋话阒傅乇硪韵?00 m內(nèi)所賦存的地?zé)豳Y源。我國地?zé)豳Y源分布廣泛、儲量豐富，淺層地?zé)衢_發(fā)前景十分廣闊［2］。就河南省而言，18個省轄市內(nèi)皆有可采潛力巨大的中低溫地?zé)豳Y源［3］，全省地?zé)峥刹闪黧w量為9.222 393×108m3／a，可利用熱能量為196.46×1012kJ／a，約合標(biāo)準(zhǔn)煤6.703 5×106t／a［4］。目前，河南省開采使用的地?zé)豳Y源主要來自被埋藏于新生界與前新生界熱儲層中的孔隙水、巖溶水和裂隙水［5］。在開發(fā)利用中，還存在勘查研究程度低與開發(fā)深度淺、地?zé)嵫h(huán)演化機(jī)理還沒有理清、地?zé)衢_發(fā)利用模式相對簡單等問題。

近年來，隨著計算機(jī)行業(yè)不斷發(fā)展，數(shù)值模擬在地下水領(lǐng)域被廣泛利用。Visual Mod flow是目前被國內(nèi)外學(xué)者普遍認(rèn)可并使用的專業(yè)化地下水流系統(tǒng)模擬軟件［6］，我國最早在20世紀(jì)90年代開始應(yīng)用［7］，在地下水資源量評價、礦井水涌水量預(yù)測、地下水位動態(tài)研究等方面得到廣泛應(yīng)用［8］。刁維杰等［9］通過建立地下水?dāng)?shù)值模型，解決了濰坊市北部地區(qū)的水資源配置問題；高學(xué)平等［10］為定量分析黎河下游的水量平衡情況，構(gòu)建了耦合HEC-RAS與Mod flow的河流-地下水系統(tǒng)模型，得到了令人滿意的結(jié)果；劉基等［11］利用Mod flow對礦井涌水量的模擬與預(yù)測進(jìn)行了深入探討；李文雅等［12］為評價海勃灣水庫的地質(zhì)環(huán)境問題，基于Visual Mod flow對蓄水后的海勃灣水庫地下滲流場進(jìn)行了模擬，取得了良好效果；魏傳云等［13］、李巍 等［14］、饒磊等［15］均將Visual Mod flow應(yīng)用于地下水污染防治，為研究地下水環(huán)境評價提供可靠依據(jù)。但在地?zé)嵫芯恐校ㄌ貏e是在地?zé)崴Y源儲量豐富的河南地區(qū)），數(shù)值模擬應(yīng)用研究的成果較少［16］，為探索新的有效途徑，本文在分析中引入此法。時間序列法［17］是利用井孔水位降深隨時間變化的數(shù)據(jù)定量計算地下水可開采量、預(yù)測地下水位動態(tài)變化的實用性理論，在地下水動力學(xué)領(lǐng)域中廣泛使用，也是分析地下淺層熱儲可開采量的普適性方法。本文利用Visual Mod flow模擬軟件建立地下淺層熱儲數(shù)值模型，并與傳統(tǒng)地下水動力學(xué)法進(jìn)行對比，深入探究兩種分法的優(yōu)劣和結(jié)合的優(yōu)勢，對河南省汝州市溫泉鎮(zhèn)淺層地?zé)崽锏目砷_采量進(jìn)行科學(xué)預(yù)測，以期為該區(qū)域的淺層地?zé)崽锖侠黹_發(fā)提供依據(jù)。

1 研究區(qū)概況及地?zé)嵫h(huán)機(jī)理

研究區(qū)位于汝州市西部，南臨汝河，北望箕山，西傍廣成湖（澗山口水庫），東為汝河沖積平原，占地面積約8.2 km2。區(qū)域內(nèi)地層自老至新出露為太古界、中元古界、寒武系、石炭系、二疊系、三疊系、新近系、第四系。溫泉鎮(zhèn)地?zé)崽锷w層由第四系松散層構(gòu)成，平均厚度為12 m。地下含水系統(tǒng)包括松散巖類孔隙水、碳酸鹽巖溶裂隙溶洞水、碎屑巖類裂隙孔隙水和巖漿巖及變質(zhì)巖裂隙水4種類型，其中碳酸鹽巖溶裂隙水是溫泉鎮(zhèn)地?zé)崴闹饕獊碓?，主要賦存于灰?guī)r巖溶裂隙含水層中，該含水層巖溶裂隙發(fā)育，富水性強(qiáng)，水量豐富。區(qū)內(nèi)熱水井水溫均大于30℃小于70℃，屬于低溫地?zé)崽?，主要?yīng)用于采暖、理療和洗浴。

溫泉鎮(zhèn)地?zé)崽镂挥诒蔽飨驍嗔褞В‵1，F(xiàn)2，F(xiàn)3）與北東向構(gòu)造（F6，F(xiàn)7）的交匯部位（圖1）。地?zé)崽锼疅峄顒訃?yán)格受斷裂構(gòu)造（斷裂）控制，具有斷裂深循環(huán)型（對流型）地?zé)嵯到y(tǒng)特征。圖2主要反映了該區(qū)域深循環(huán)地?zé)嵯到y(tǒng)熱量和質(zhì)量傳遞過程。從地下較遠(yuǎn)處補(bǔ)給深部熱水（Tr）沿通道（深度H）上升，上升熱水在一定深度（Hi）和來自透水層中冷水（T0）相遇混合，混合水再沿通道上升（溫度并不斷散失），溢出地表便形成溫泉（Ts）。從區(qū)域地質(zhì)構(gòu)造條件分析，北西向斷裂帶為九皋山-溫泉街?jǐn)嗔押托掳?平頂山斷裂在此交匯復(fù)合形成，九皋山-溫泉街?jǐn)嗔押托掳?平頂山均為區(qū)域性的控?zé)釘嗔?。東西與北西向（F1F3）主干斷裂以導(dǎo)水為主、導(dǎo)熱為輔的復(fù)合斷裂帶，主要將西北部碳酸鹽巖水輸送到溫泉鎮(zhèn)熱儲層內(nèi)。北東向（F6）以導(dǎo)熱為主、導(dǎo)水為輔的次級斷裂帶，控制著溫泉的形成、儲存與出露。

圖1 研究區(qū)地?zé)崽锂惓^(qū)分布范圍Fig.1 Distribution range of geothermal field anomalies in the study area

圖2 研究區(qū)深循環(huán)水熱系統(tǒng)示意圖Fig.2 Schematic diagram of deep-circulation hydrothermal system in the study area

2 地?zé)峥砷_采量預(yù)測

2.1 時間序列法預(yù)測

2.1.1 時間序列法

依據(jù)《地?zé)豳Y源地質(zhì)勘查規(guī)范》（GB11615-2010）規(guī)定，單個地?zé)峋?，確定開采量使用的壓力降低值一般不大于0.3 MPa（約合30.6 m）。而對于多井開采的地?zé)崽?，主要以地?zé)崽飪?nèi)代表性監(jiān)測井保持一定水位年降速條件下的地?zé)崽镩_采量作為一定時限內(nèi)的可開采量。首先，通過群孔抽水試驗資料，對短期開采狀態(tài)下的主要井孔水位進(jìn)行模擬，得出基本適應(yīng)于本地?zé)崽锏乃牡刭|(zhì)參數(shù)；其次，在一定時期內(nèi)推算出在適當(dāng)降深情況下的抽水量，以此抽水量作為可采地?zé)豳Y源。

2.1.2 預(yù)測步驟

（1）數(shù)學(xué)模型及邊界條件概化。根據(jù)研究區(qū)物探解譯成果，溫泉鎮(zhèn)地?zé)豳x存區(qū)是沿北西向斷裂帶（F1—F3）為南部北部邊界，F(xiàn)6斷層以西區(qū)域，F(xiàn)6斷層為隔水邊界，地下水位埋深4.7～7.7 m。其余地段地下水儲水能力差。為便于計算，以F6斷層為界，數(shù)學(xué)模型概化為半無限承壓含水層。據(jù)此，采用干擾井群非穩(wěn)定流公式計算允許開采量及降深。

（2）計算公式的選擇。根據(jù)概化的邊界條件，按映射原理以F6直線隔水邊界為對稱軸進(jìn)行映射，利用承壓完整井干擾非穩(wěn)定流計算，其計算公式為

式中：s為任意點處的降深，m；Qi為第i個井的出水量，m3／d，包括虛擬井（水量為正）；T為導(dǎo)水系數(shù)，m2／d；W（ui）為井函數(shù)；u′i為虛擬井自變量；ui為井函數(shù)自變量，ri為計算點到第i個井（含虛擬井）的距離，m；S為彈性釋水系數(shù)；t為抽水延續(xù)時間，d。

（3）布井方案的確定。根據(jù)溫泉鎮(zhèn)地?zé)岣凰畮У奈恢?，結(jié)合研究區(qū)水文地質(zhì)條件，選擇F2以南、F3以北、F6以西、F7以東區(qū)域布井（圖3），單井出水量分別為70，43，30m3／h，總開采量3432m3／d。虛擬井3口，水量3 432 m3／d。

（4）參數(shù)驗證。將前期推算出的研究區(qū)水文地質(zhì)參數(shù)代入干擾井群公式反演，以2018年9—10月群孔抽水試驗資料為依據(jù)，抽水量分別為70，43，30 m3／h，虛擬井水量與抽水井相對應(yīng)，以擬合14 d ZK1孔抽水時降深為主。經(jīng)模擬反演，當(dāng)T＝253.07 m2／d，彈性釋水系數(shù)S＝0.006 92時與抽水資料吻合（計算降深15.81 m，抽水試驗降深16.508 m，扣除井損0.707 m，實際降深為15.801 m）。

（5）最大降深與可允許開采量的計算。考慮溫泉鎮(zhèn)地?zé)釋儆趲顭醿?，地?zé)豳Y源有限，將反演參數(shù)代入干擾井群公式計算，初步按100 d計算，累計降深30 m左右。經(jīng)計算，開采100 d后，ZK1井為水位下降中心區(qū)，水位降深33.395 m（表1）。

由表1可知，為滿足長期開采的需要，當(dāng)群井方案中的抽水量為3 432 m3／d且以100 d為期進(jìn)行計算時，ZK1井的降深基本穩(wěn)定不變。故而可得溫泉鎮(zhèn)地?zé)崽锏目稍试S開采量為3 432 m3／d。

2.2 基于Mod flow的地?zé)崃黧w數(shù)值模擬

2.2.1 溫泉鎮(zhèn)地?zé)崽锼牡刭|(zhì)概念模型

模型構(gòu)建范圍為汝州市溫泉鎮(zhèn)淺低溫地?zé)崽?。地?zé)崽飮?yán)格受北西向斷裂（F1—F5）和北東向斷裂（F6—F7）控制，其中F2和F3穿越了寒武系灰?guī)r，構(gòu)成地?zé)崴幕顒訄鏊?，因此沿北西向斷裂帶（F2—F3）和北東向構(gòu)造（F6—F7）之間形成了地?zé)岙惓^(qū)（圖1）。為準(zhǔn)確模擬溫泉鎮(zhèn)地?zé)崽锏牡叵滤疇顟B(tài)，模擬過程將南北部邊界分別向外延伸至F1和F3。因此，模型的南部邊界為F1，北部邊界為F3，東部邊界為F6，西部邊界為F7。

模型主要以賦存在灰?guī)r巖溶裂隙含水層中的碳酸鹽巖溶裂隙水進(jìn)行構(gòu)建。該含水層受構(gòu)造控制明顯，處在斷層帶上，多個層段巖石破碎或裂隙較發(fā)育，水文地質(zhì)條件較好，水文地質(zhì)參數(shù)隨空間變化明顯，垂直方向與水平方向滲透系數(shù)不同，故概化為非均質(zhì)、各向異性。含水層上覆第四系沖洪積物，并有砂質(zhì)黏土巖和黏土層，厚度較薄，為7～13 m，構(gòu)成相對隔水層，地下水處于承壓狀態(tài)。綜上所述，將其概化為非均質(zhì)各向異性的承壓非穩(wěn)定流模型。

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

根據(jù)上述的水文地質(zhì)概念模型，可用以下數(shù)學(xué)模型進(jìn)行描述：

式中：kxx，kyy，kzz為沿x，y，z坐標(biāo)軸方向的滲透系數(shù)，m·d－1；h為點（x，y，z）在t時刻的水位標(biāo)高，m；W為源匯項，L·d－1；Ss為單位儲水系數(shù)，L·m－1；t為時間，d；Ω為承壓區(qū)域；S1為第一類邊界；S2為第二類邊界；φ（x，y，z，t）為第一類邊界上水位標(biāo)高；q（x，y，z，t）為流量在時間和空間上的變化函數(shù)；n為二類邊界內(nèi)法線向量。

2.2.3 研究區(qū)空間離散

應(yīng)用地下水?dāng)?shù)值模擬軟件Visual Mod flow對研究區(qū)進(jìn)行等距或不等距長方體網(wǎng)格剖分，得到空間離散單元，為提高模擬的精度，對觀測井、抽水井、邊界等區(qū)域進(jìn)行局部加密，最終在平面上剖分成92×92的矩形網(wǎng)格單元，垂向上為一層，即碳酸鹽巖溶裂隙含水層，共計8 464個網(wǎng)格單元，其中有效單元格3 182個，無效單元格5 282個。

2.2.4 邊界、源匯項處理

南部邊界：F1斷裂帶位于溫泉鎮(zhèn)南200 m，該斷層將寒武系灰?guī)r阻斷，缺少熱活動必備的空間條件，起隔水作用，故設(shè)為阻水邊界；北部邊界：F3斷裂帶同樣將寒武系灰?guī)r阻斷，起隔水作用，設(shè)為阻水邊界；東部邊界：F6斷層切斷了所有北西向斷裂，使寒武系在斷裂以東消失，阻斷了熱水沿F3斷裂向東南的通道，設(shè)為阻水邊界；西部邊界：F7斷裂帶接受地下水補(bǔ)給，設(shè)為一般水位邊界（GHB）［18］。

本區(qū)地下水流的補(bǔ)給來源為斷裂帶降水入滲補(bǔ)給和側(cè)向徑流補(bǔ)給，以大氣降水的垂向補(bǔ)給為主。大氣降水經(jīng)灰?guī)r巖溶裂隙、斷破碎帶滲入地下深部，匯入F2，F(xiàn)3斷裂帶，入滲系數(shù)取0.5。通過Visual Mod flow中的RECHARGE程序包，可計算降水入滲補(bǔ)給量。側(cè)向徑流補(bǔ)給量主要在西部邊界，將其概化為補(bǔ)給邊界。地下水排泄以人工鉆井排泄為主要途徑，排泄項采用抽水井的方式表示。

2.2.5 水文地質(zhì)參數(shù)分區(qū)

結(jié)合溫泉鎮(zhèn)的水文地質(zhì)資料及前人研究成果，將研究區(qū)劃分為5個參數(shù)分區(qū)（圖4），并賦予儲水系數(shù)Ss和滲透系數(shù)k初值，作為模型計算的基礎(chǔ)。

圖4 參數(shù)分區(qū)圖Fig.4 Parameter partition map

2.2.6 模型的識別與驗證

根據(jù)群孔抽水試驗資料，選取2018年9月2日至2018年9月16日作為模型的校正、識別階段，整個模擬期共分3個制度期，每個制度期又分10個計算時間步長。此次群孔抽水試驗，共設(shè)有抽水井5個，分別為ZK1號、3號（中井）、4號（東井）、療2和療3。

通過運行該模型，根據(jù)觀測井預(yù)測水位值與實測值的擬合情況，不斷地試算參數(shù)，直至兩者擬合較好，得到各分區(qū)的參數(shù)值。經(jīng)校正、識別，各分區(qū)的參數(shù)識別結(jié)果參見表2。

表2 各參數(shù)分區(qū)參數(shù)識別結(jié)果Tab.2 Parameter identification results of each parameter partition

各水位觀測井?dāng)M合精度如圖5所示。

由圖5可以看到，在抽水試驗過程中，各水位觀測孔均分布在95%的置信區(qū)間內(nèi)，達(dá)到了模型識別和驗證的要求。模擬結(jié)果表明，水文地質(zhì)概念模型基本合理，試算結(jié)束后的水文地質(zhì)參數(shù)和各源匯項與實際情況基本相符，可以認(rèn)為本次所建立的數(shù)值模型基本反映了模擬區(qū)的地?zé)崃黧w運動規(guī)律?？蓪⒅畱?yīng)用于地?zé)崃黧w流場和可開采量動態(tài)預(yù)測等相關(guān)研究。

圖5 地下水位觀測井?dāng)M合精度Fig.5 Groundwater level fitting accuracy of the observation well

2.2.7 研究區(qū)可開采量預(yù)測

采用數(shù)值模擬方法確定地下水開采量，應(yīng)在地下水總補(bǔ)給和總排泄量均衡的前提下，以使地下水得到充分的開采利用，水位不出現(xiàn)持續(xù)下降為標(biāo)準(zhǔn)。通過調(diào)節(jié)抽水井抽水強(qiáng)度，使水位只在某一范圍變動或趨于穩(wěn)定，沒有持續(xù)下降的趨勢，將這種條件下的開采量作為最大可開采量。綜上所述，在長期開采條件下，將各水位觀測井水位逐漸趨于穩(wěn)定時的模擬開采量作為溫泉鎮(zhèn)地?zé)崽锏淖畲罂砷_采量。以此為約束條件，對研究區(qū)未來5年（2018年9月17日至2023年9月17日）的地下水最大可開采量和地下水位動態(tài)變化進(jìn)行模擬預(yù)測（圖6）。

圖6 模擬期各水位觀測井水位變化曲線Fig.6 Water level change curves of each water level observation well in the simulation period

模擬結(jié)果表明，當(dāng)ZK1井開采量為1 200 m3／d，1號井開采量為480 m3／d，3號井開采量為720 m3／d，總開采量為2 400 m3／d時，各觀測井水位變化逐漸趨于穩(wěn)定，各水位觀測井水位變化情況如圖6所示。因此，確定溫泉鎮(zhèn)地?zé)崽锏責(zé)崃黧w的最大可開采量約為2 400 m3／d。

2.3 討論

賦存在地下淺層熱儲中的地?zé)崃黧w受斷層邊界、源匯項等復(fù)雜地質(zhì)環(huán)境的影響，其運移循環(huán)機(jī)理是多種復(fù)雜變量相互疊加作用后的綜合顯現(xiàn)，具有特定的水文地質(zhì)物理機(jī)制。時間序列法是Thesis公式在地下水問題研究中的經(jīng)典應(yīng)用，也可作為研究淺層地?zé)崴畣栴}的一種方法。它將抽象復(fù)雜的動態(tài)物理過程轉(zhuǎn)變成具有概化性的數(shù)學(xué)模型，通過簡潔的數(shù)學(xué)方法對地?zé)衢_采后的水位作出評價；Visual Mod flow數(shù)值模擬法就是在嚴(yán)格的物理成因基礎(chǔ)上，將研究區(qū)的各類限制條件輸入到所構(gòu)建的模型中，以此對地下熱儲層進(jìn)行研究分析。前者利用計算機(jī)模擬地?zé)釀討B(tài)變化特征，更加生動、詳實，且以模型精度為保證，具有更高的可信度。但需要大量而準(zhǔn)確的野外勘查數(shù)據(jù)為支撐，適用于后期室內(nèi)階段的分析研究；后者將復(fù)雜的地下水流問題轉(zhuǎn)變?yōu)閿?shù)學(xué)建模理論，僅從井孔水位降深隨時間變化關(guān)系的角度出發(fā)進(jìn)行探究，減少了對水文地質(zhì)資料的依賴程度，為預(yù)測水位變化提供便捷、有效的途徑。但真實的地?zé)崃黧w處于受各類外界變量影響的非線性地質(zhì)環(huán)境系統(tǒng)中，此種方法考慮因素單一，不能完全反映地?zé)崴鞯倪\動狀態(tài)，只能對可開采量進(jìn)行簡單粗略地估計，更適合開采現(xiàn)場的試驗觀測。本文將室內(nèi)解析反演法與室內(nèi)數(shù)值模擬法相結(jié)合，通過兩種方法的互補(bǔ)優(yōu)勢，為地?zé)豳Y源的實際生產(chǎn)開發(fā)提供有效的參考依據(jù)。

近年來，隨著計算機(jī)信息技術(shù)在地?zé)豳Y源領(lǐng)域中的不斷應(yīng)用，相比于利用傳統(tǒng)時間-降深序列的可開采量預(yù)測，涵蓋多角度地質(zhì)物理特性的數(shù)值模擬法對淺層熱儲流體的模擬將會得到長足發(fā)展。

3 結(jié)論

（1）汝州市溫泉鎮(zhèn)地?zé)崽餅闇\層低溫地?zé)崽?，地?zé)犷愋蛯贁嗔焉钛h(huán)型地?zé)嵯到y(tǒng)。熱儲層由石炭系、寒武系灰?guī)r破碎帶組成，地下熱水為裂隙承壓水。

（2）采用時間序列法得到研究區(qū)的可允許開采量為3 432 m3／d；利用Visual Mod flow數(shù)值模擬法預(yù)測未來5年水位降深達(dá)到穩(wěn)定狀態(tài)后，地?zé)崃黧w的最大可開采量為2 400 m3／d。結(jié)合研究區(qū)水文地質(zhì)條件，汝州市溫泉鎮(zhèn)淺層地?zé)崽锏牡責(zé)崃黧w的最大可開采量應(yīng)處于2 400～3 432 m3／d之間。淺層地?zé)嶂袘?yīng)用Visual Mod flow數(shù)值模擬法，實現(xiàn)了現(xiàn)場試驗技術(shù)與室內(nèi)研究方法的有機(jī)結(jié)合，可為地?zé)峥砷_采量的確定提供有效支撐。