基于Mod flow和時間序列法的汝州市溫泉鎮(zhèn)地熱可開采量預測

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

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

0 引言

地熱資源是一種清潔、低碳的可再生能源。地熱資源的開發(fā)利用可有效優(yōu)化能源供給結構，對解決全球化的生態(tài)環(huán)境問題具有深遠影響。按照埋藏深度，可將地熱資源劃分為淺層地熱、水熱型地熱和干熱巖［1］。淺層地熱一般指地表以下200 m內所賦存的地熱資源。我國地熱資源分布廣泛、儲量豐富，淺層地熱開發(fā)前景十分廣闊［2］。就河南省而言，18個省轄市內皆有可采潛力巨大的中低溫地熱資源［3］，全省地熱可采流體量為9.222 393×108m3／a，可利用熱能量為196.46×1012kJ／a，約合標準煤6.703 5×106t／a［4］。目前，河南省開采使用的地熱資源主要來自被埋藏于新生界與前新生界熱儲層中的孔隙水、巖溶水和裂隙水［5］。在開發(fā)利用中，還存在勘查研究程度低與開發(fā)深度淺、地熱循環(huán)演化機理還沒有理清、地熱開發(fā)利用模式相對簡單等問題。

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

1 研究區(qū)概況及地熱循環(huán)機理

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

溫泉鎮(zhèn)地熱田位于北西向斷裂帶（F1，F(xiàn)2，F(xiàn)3）與北東向構造（F6，F(xiàn)7）的交匯部位（圖1）。地熱田水熱活動嚴格受斷裂構造（斷裂）控制，具有斷裂深循環(huán)型（對流型）地熱系統(tǒng)特征。圖2主要反映了該區(qū)域深循環(huán)地熱系統(tǒng)熱量和質量傳遞過程。從地下較遠處補給深部熱水（Tr）沿通道（深度H）上升，上升熱水在一定深度（Hi）和來自透水層中冷水（T0）相遇混合，混合水再沿通道上升（溫度并不斷散失），溢出地表便形成溫泉（Ts）。從區(qū)域地質構造條件分析，北西向斷裂帶為九皋山-溫泉街斷裂和新安-平頂山斷裂在此交匯復合形成，九皋山-溫泉街斷裂和新安-平頂山均為區(qū)域性的控熱斷裂。東西與北西向（F1F3）主干斷裂以導水為主、導熱為輔的復合斷裂帶，主要將西北部碳酸鹽巖水輸送到溫泉鎮(zhèn)熱儲層內。北東向（F6）以導熱為主、導水為輔的次級斷裂帶，控制著溫泉的形成、儲存與出露。

圖1 研究區(qū)地熱田異常區(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 地熱可開采量預測

2.1 時間序列法預測

2.1.1 時間序列法

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

2.1.2 預測步驟

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

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

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

（3）布井方案的確定。根據溫泉鎮(zhèn)地熱富水帶的位置，結合研究區(qū)水文地質條件，選擇F2以南、F3以北、F6以西、F7以東區(qū)域布井（圖3），單井出水量分別為70，43，30m3／h，總開采量3432m3／d。虛擬井3口，水量3 432 m3／d。

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

（5）最大降深與可允許開采量的計算?？紤]溫泉鎮(zhèn)地熱屬于帶狀熱儲，地熱資源有限，將反演參數代入干擾井群公式計算，初步按100 d計算，累計降深30 m左右。經計算，開采100 d后，ZK1井為水位下降中心區(qū)，水位降深33.395 m（表1）。

由表1可知，為滿足長期開采的需要，當群井方案中的抽水量為3 432 m3／d且以100 d為期進行計算時，ZK1井的降深基本穩(wěn)定不變。故而可得溫泉鎮(zhèn)地熱田的可允許開采量為3 432 m3／d。

2.2 基于Mod flow的地熱流體數值模擬

2.2.1 溫泉鎮(zhèn)地熱田水文地質概念模型

模型構建范圍為汝州市溫泉鎮(zhèn)淺低溫地熱田。地熱田嚴格受北西向斷裂（F1—F5）和北東向斷裂（F6—F7）控制，其中F2和F3穿越了寒武系灰?guī)r，構成地熱水的活動場所，因此沿北西向斷裂帶（F2—F3）和北東向構造（F6—F7）之間形成了地熱異常區(qū)（圖1）。為準確模擬溫泉鎮(zhèn)地熱田的地下水狀態(tài)，模擬過程將南北部邊界分別向外延伸至F1和F3。因此，模型的南部邊界為F1，北部邊界為F3，東部邊界為F6，西部邊界為F7。

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

2.2.2 數學模型

根據上述的水文地質概念模型，可用以下數學模型進行描述：

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

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

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

2.2.4 邊界、源匯項處理

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

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

2.2.5 水文地質參數分區(qū)

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

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

2.2.6 模型的識別與驗證

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

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

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

各水位觀測井擬合精度如圖5所示。

由圖5可以看到，在抽水試驗過程中，各水位觀測孔均分布在95%的置信區(qū)間內，達到了模型識別和驗證的要求。模擬結果表明，水文地質概念模型基本合理，試算結束后的水文地質參數和各源匯項與實際情況基本相符，可以認為本次所建立的數值模型基本反映了模擬區(qū)的地熱流體運動規(guī)律。可將之應用于地熱流體流場和可開采量動態(tài)預測等相關研究。

圖5 地下水位觀測井擬合精度Fig.5 Groundwater level fitting accuracy of the observation well

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

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

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

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

2.3 討論

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

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

3 結論

（1）汝州市溫泉鎮(zhèn)地熱田為淺層低溫地熱田，地熱類型屬斷裂深循環(huán)型地熱系統(tǒng)。熱儲層由石炭系、寒武系灰?guī)r破碎帶組成，地下熱水為裂隙承壓水。

（2）采用時間序列法得到研究區(qū)的可允許開采量為3 432 m3／d；利用Visual Mod flow數值模擬法預測未來5年水位降深達到穩(wěn)定狀態(tài)后，地熱流體的最大可開采量為2 400 m3／d。結合研究區(qū)水文地質條件，汝州市溫泉鎮(zhèn)淺層地熱田的地熱流體的最大可開采量應處于2 400～3 432 m3／d之間。淺層地熱中應用Visual Mod flow數值模擬法，實現(xiàn)了現(xiàn)場試驗技術與室內研究方法的有機結合，可為地熱可開采量的確定提供有效支撐。