宋前進，王 剛，許一川，程 磊

(1.中國地質大學(北京) 工程技術學院，北京 100083；2.河南省地質礦產(chǎn)勘查開發(fā)局第二地質環(huán)境調查院，河南 鄭州 450053；3.河南財經(jīng)政法大學 工程管理與房地產(chǎn)學院，河南 鄭州 450046)

地熱資源是指貯存在地球內部的可再生熱能，一般集中分布在構造板塊邊緣一帶，起源于地球的熔融巖漿和放射性物質的衰變。根據(jù)地熱資源所處深度，可將其劃分為淺層地溫能、深層地熱能和干熱巖等類型[1]。深層地熱能需有特定的儲藏條件，即熱源、導熱斷裂、控熱斷裂、有效熱儲層和蓋層等要素，地熱水流場也與這些要素有關，其中熱儲層特征對地熱水流場的影響不容忽視[2-5]。文慧儉等[6]認為有效熱儲層是具有一定埋深和厚度、可使地層水采出地面時仍有40～60℃的具有商用價值的含水巖層；李雪英等[7]對松遼盆地熱儲層滲流模型開展研究，基于質量和能量守恒原理、達西滲流規(guī)律，建立熱儲層流場的滲流數(shù)學模型；張萌等[8]研究浙江省地熱資源的賦存特征，詳細分析熱儲層特征；李海泉等[9]從區(qū)域地質構造角度研究臨汾盆地曲沃地熱田熱儲特征，明確地熱田的熱儲邊界；馬月花等[10]以青海共和盆地地熱資源勘探開發(fā)為例，從地層巖性、地質構造角度分析熱儲層特征對地熱水流場的重要作用。

深層地熱能是通過開發(fā)其載體得以利用，開采時會改變地熱水的補徑排條件，進而改變熱儲層地熱水的流場狀態(tài)。長期不合理開采會造成地熱水水位持續(xù)下降，形成一定規(guī)模的降落漏斗，進而引發(fā)區(qū)域性地面沉降等環(huán)境地質問題。筆者以蘭考縣深層地熱能開發(fā)為例，分析了區(qū)域地質條件和深層地熱的儲藏情況，并分區(qū)論述地熱資源量及開發(fā)利用前景，運用數(shù)值模擬方法預測深層地熱能開發(fā)對地熱水流場造成的影響，以期為蘭考縣和周邊地區(qū)地熱資源開發(fā)及環(huán)境保護提供參考。

1 區(qū)域地質與地熱開采狀況

1.1 地質與構造

蘭考縣位于河南省東部，地熱清潔能源開發(fā)利用在全省范圍內起到了示范性作用。該縣地處黃淮平原腹地，屬流水堆積地貌，新生代以來，區(qū)內沉積巨厚新生代地層，包括古近系、新近系和第四系。古近系分為沙河街組和東營組，頂板埋深約2 000 m，底板埋深4 000～6 500 m，巖性有泥巖、砂巖、含礫砂巖及角礫狀灰?guī)r；新近系分為館陶組和明化鎮(zhèn)組；館陶組(N1)為泥巖、砂質泥巖、塊狀砂巖，平均厚度650 m；明化鎮(zhèn)組(N2)為砂質泥巖與粉細砂巖互層，底部為含礫砂巖或砂礫巖，平均厚度1 100 m；第四紀地層底板埋深300～360 m，巖性為粉細砂、粉質黏土及黏土。

根據(jù)地質構造區(qū)劃，蘭考縣屬華北地臺區(qū)，位于秦嶺東西向構造體系東端和新華夏構造體系的復合部位。該區(qū)被新鄉(xiāng)-商丘大斷裂(F3)分為南北兩部分，北部在東濮凹陷內，南部在開封凹陷內，基底被NE、NW和EW向3組斷裂切割為不同的次級斷塊構造，如圖1所示。構造活動總體特征為斷陷，不同斷塊的升降運動差異性明顯。

1.2 地熱開采歷史、現(xiàn)狀及存在問題

自20世紀80年代，蘭考縣開始零星開發(fā)地熱，在2000年前后形成規(guī)模。2013年前主要是利用第四系(Q)和新近系明化鎮(zhèn)組(N2)的地熱能，地熱井深度1 000～1 300 m，水溫40℃左右。2013年，河南省中石化新星石油有限公司在蘭考油田六社區(qū)內施工1口探采結合試驗井，成井深度1 980 m，單井出水量120 m3/h，出孔水溫72℃，開采對象為新近系館陶組熱儲層水。隨后，蘭考油田六社區(qū)相繼又施工3口井深2 000 m的地熱井。截至2019年12月底，區(qū)內井深1 942～2 270 m的地熱井超過24口，深度小于1 300 m的地熱井逐步被停用。新近系館陶組(N1)熱儲成為區(qū)內地熱能的主要開發(fā)對象。

熱儲層地熱水的天然流場受區(qū)域地質構造、地層條件控制[11-13]。隨著地熱開發(fā)，人工開采成為熱儲層熱水的主要排泄方式，對流場的影響也越來越大。前期開采第四系和新近系明化鎮(zhèn)組地熱能時，熱水用完后基本排入河流等地表水體，造成很大資源浪費。地熱開發(fā)到一定規(guī)模后，熱儲層熱水的排泄量遠大于補給量，造成地下水水位持續(xù)下降。根據(jù)蘭考縣三義寨鄉(xiāng)三義寨村地熱井的水位觀測資料，新近系明化鎮(zhèn)組熱儲層水位下降情況如圖2所示，下降速度驚人。水位持續(xù)下降會形成不同規(guī)模的降落漏斗，導致原有的流場遭到破壞。同時，熱儲層水位降低改變了熱儲層及上覆地層的應力狀況，進而引發(fā)地面沉降等環(huán)境問題。

圖1 區(qū)域地質構造Fig.1 Regional geological structure

圖2 新近系明化鎮(zhèn)組熱儲層水位下降曲線Fig.2 Decline of heat storage head at Neogene Minghua Town Formation

近幾年，隨著新近系館陶組地熱能的開發(fā)利用，第四系和新近系明化鎮(zhèn)組地熱的開采井逐漸被停用。2018年后，區(qū)內地熱開發(fā)利用方式也發(fā)生變化，開始采用封閉循環(huán)利用熱能，將取熱后回水同層回灌熱儲層。根據(jù)應用情況，一對一抽采回灌時回灌率普遍較低，一般在40%～50%。

2 區(qū)域地熱資源特征

2.1 區(qū)域熱儲層

蘭考縣深層地熱能屬沉積盆地熱儲類型[14]，熱源供給為大地熱流傳導，即地球深部的熱能傳輸至地表。區(qū)內新生界地層厚度大，具備形成熱儲層的條件，淺部第四系發(fā)育有多層黏土和砂質黏土層，具有保溫隔熱作用，有利于地熱的富集與儲存。區(qū)內主要有第四系(Q)、新近系明化鎮(zhèn)組(N2)、新近系館陶組(N1)及古近系等熱儲層。第四系(Q)熱儲層深度300～350 m，水溫略高于25℃；明化鎮(zhèn)組(N2)熱儲層深度350～1 300 m，水溫30～55℃；館陶組(N1)熱儲層深度1 500～2 200 m，水溫80～90℃。目前，本區(qū)古近系及其以下的地熱資源尚未開發(fā)。

2.2 構造影響

區(qū)內發(fā)育有聊城-蘭考斷裂(F1)、鄭州-蘭考斷裂(F2)、新鄉(xiāng)-商丘斷裂(F3)等多條區(qū)域性斷裂，切穿地殼，活動強烈[15]。這些斷裂形成良好的地下熱水運移通道和儲存空間，使深部熱能可通過F3、F7等斷裂傳至新近系館陶組(N1)、明化鎮(zhèn)組(N2)及第四系的熱儲層中，如圖3所示。傳導至N1、N2及Q的熱能被上覆第四系黏土層、砂質黏土層等阻隔得以大量儲存，形成可開采的深層地熱資源。熱儲層的地熱資源量由地層結構和沉積特征決定，與導熱斷裂發(fā)育情況密切相關。根據(jù)地熱調查，區(qū)內淺層地下水溫度分布的特點為基底斷裂構造越發(fā)育，溫度越高、地熱異常越明顯。

2.3 熱儲層的差異及影響因素

圖3 蘭考縣熱儲層結構Fig.3 Thermal reservoir structure of Lankao County

蘭考縣深層地熱水的補給源主要為開封凹陷西部山區(qū)大氣降水，地下水排泄方式以側向徑流和人工開采為主。研究區(qū)深層地熱能在地表無地熱顯示，自恒溫帶以下溫度隨深度的增加而升高。根據(jù)蘭考職業(yè)技術學院地熱井施工時的測井資料，地層溫度隨深度的增加曲線如圖4所示。研究區(qū)目前所開采的深層地熱能有第四系(Q)熱儲層、明化鎮(zhèn)組(N2)熱儲層和館陶組(N1)熱儲層。第四系熱儲層深度300～350 m，新近系明化鎮(zhèn)組熱儲埋深1 000～ 1 300 m，新近系館陶組熱儲埋深1 900～2 400 m，3個熱儲層的地下水在地質、地溫、水文地質、水化學特征等方面明顯不同，為3個相互獨立的熱儲層。由于各熱儲層之間存在較厚的隔水層，各層之間的垂向水力聯(lián)系不明顯，不同熱儲層的水位高度差異較大。例如蘭考縣三義寨鄉(xiāng)三義寨村礦泉水井，成井深度650 m，開采層埋深463.6～644.35 m，水位埋深為70.12 m；蘭考縣天泉洗浴廣場的地熱井，成井深度1 161.4 m，開采層埋深1 087.4～1 153.3 m，水位埋深27 m；蘭考縣油田六社區(qū)1號地熱井，成井深度1 980 m，開采層埋深1 511.4～1 972.3 m，水位埋深為60.69 m。

圖4 蘭考縣地層溫度隨深度變化曲線Fig.4 Temperature-increasing curve of strata with depth in Lankao County

研究區(qū)內各熱儲層埋藏條件不同，熱儲量和資源量差別較大，不同構造部位熱儲層的水文地質參數(shù)也有所不同，見表1。處在不同的構造部位，受到斷裂構造的影響，導致這些地帶熱儲層裂隙較發(fā)育，形成地熱水良好的運移通道和儲存空間，熱儲層的富水性與斷裂構造的關系極為密切。油田六社區(qū)地熱井在一抽一回狀態(tài)下，其回灌量可達114.8 m3/h，回灌系統(tǒng)穩(wěn)定壓力水頭高19.65 m，回灌效果比較好，推測其原因也與斷裂附近的裂隙發(fā)育密度密切相關。

另外，熱儲層同層之間地下水的流動性提高了巖體導熱效果，使地熱梯度值變小，地溫場的垂向變化一定程度上也受地層巖性和地下水活動情況的影響。

2.4 地熱分區(qū)及地熱資源量

針對目前研究區(qū)所開發(fā)的3個深層熱儲層，采用GDP-32Ⅱ多功能電法儀對全區(qū)280 km2進行可控源音頻大地電磁測深[16]，查明區(qū)內的導熱斷裂分布、熱儲層發(fā)育、隔水邊界等條件，并通過地質勘探方法進行印證。區(qū)內不同部位熱儲層的厚度、埋深及地熱資源量有所差異。為便于計算資源量，根據(jù)熱儲層特征結合構造將全區(qū)分為4個地熱分區(qū)(Ⅰ、Ⅱ、Ⅲ、Ⅳ區(qū))，如圖5所示。

Ⅰ區(qū)面積107.78 km2，熱儲層埋深1 300～1 940 m，累計厚度201～361 m；

Ⅱ區(qū)面積57.49 km2，熱儲層埋深986～1 880 m，累計厚度50～388 m，區(qū)內新生界厚度由東南向西北逐漸減小；

表1 不同構造部位熱儲層水文地質參數(shù)Table 1 Index of thermal reservoir in different tectonic sites

圖5 蘭考縣地熱分區(qū)Fig.5 Geothermal zoning in Lankao County

Ⅲ區(qū)面積54.80 km2，熱儲層埋深1 078～1 896 m，累計厚度127～375 m，區(qū)內新生界厚度由西南向東北逐漸減?。?/p>

Ⅳ區(qū)面積59.93 km2，熱儲層埋深1 300～1 940 m，累計厚度202～357 m。

利用“熱儲法”對各區(qū)的地熱資源儲量進行計算，全區(qū)熱儲層存熱量為2.16×1019J。其中第四系總存熱(Q)為2.21×1018J，水中儲存熱量(QW)為 1.06×1018J；明化鎮(zhèn)組(N2)熱儲層存熱量(Qz)為2.92×1018J，水中儲存熱量(QW) 1.33×1018J；館陶組(N1)熱儲層存熱量(Qz)為1.65×1019J，水中存熱量(QW)為1.04×1019J，Ⅰ—Ⅳ區(qū)各層地熱資源情況見表2，包括巖體儲熱量QR。根據(jù)計算，回灌條件下館陶組地熱流體可采量為90 067 134.31 m3/a，明化鎮(zhèn)組地熱流體的可開采量約為16 468 660.80 m3/a，館陶組地熱可開采量是明化鎮(zhèn)組可采量的5.47倍。

表2 Ⅰ—Ⅳ區(qū)地熱資源情況Table 2 Ⅰ-Ⅳ region geothermal resources

2.5 針對性策略

隨著蘭考縣城市的快速發(fā)展，僅館陶組熱儲層的資源量已無法滿足經(jīng)濟社會對地熱的需求，區(qū)內古近系東營組及更深的地熱資源也應開發(fā)利用。綜合利用第四系、明化鎮(zhèn)組、館陶組及古近系東營組等多個熱儲層的資源，且合理規(guī)劃地熱開采井的平面布局[17]，是本區(qū)地熱科學、合理的開發(fā)方式。

回灌率低是國內地熱開發(fā)中的主要技術難題之一，壓力回灌是解決此問題的主要思路。為研究回灌情況，在蘭考油田六社區(qū)地熱井分別進行了小灌量加壓回灌試驗和大灌量加壓回灌試驗。小灌量回灌試驗曲線如圖6所示，回灌量78.5 m3/h時回灌系統(tǒng)穩(wěn)定壓力水頭為17.58 m；大灌量回灌試驗曲線如圖7所示，回灌量113.8 m3/h時回灌系統(tǒng)穩(wěn)定壓力水頭為19.65 m；回灌量加大44.97%，而回灌壓力水頭僅升高11.77%，即回灌系統(tǒng)加壓后，回灌率可得到有效提高。

圖6 小灌量回灌試驗曲線Fig.6 Experimental curve of small irrigation volume

圖7 大灌量回灌試驗曲線Fig.7 Experimental curve of big irrigation volume

3 熱儲層流場的數(shù)值模擬

3.1 數(shù)學模型建立

為有效控制深層地熱資源開發(fā)對區(qū)域地質環(huán)境造成的不良影響[18]，需預測地熱開采造成的熱儲層熱水流場的變化情況。根據(jù)研究區(qū)熱儲層水文地質特征及相關參數(shù)，建立數(shù)學模型模擬地熱資源開采一段時間后地熱水流場的變化情況。由于近期本區(qū)無開發(fā)古近系熱儲的規(guī)劃[14]，本次數(shù)值模擬不涉及古近系地下水流場?？紤]到區(qū)內原生地下水流場穩(wěn)定，同一含水層組地下水間存在較強的水力聯(lián)系，故將區(qū)內熱儲層地熱水流系統(tǒng)概化為均質三維穩(wěn)定地下水流系統(tǒng)，在不考慮流體密度變化的情況下，采用數(shù)學模型為：

式中：H是地下水位，m；Kx、Ky、Kz依次為x、y、z方向的滲透系數(shù)，m/d；K為邊界法線方向的滲透系數(shù)；Ss為含水層儲水率，1/m；Γ為模擬區(qū)域第一類邊界；H0是含水層初始水位，m；q為含水層第二類邊界單位面積過水斷面補給流量，m2/d；ε為源匯項強度(包括開采強度等)，d-1；Ω為滲流區(qū)域；n為滲流區(qū)邊界的單位外法線方向。

本次使用GMS(Groundwater Modeling System)10.4.5_64數(shù)值模擬軟件，水流模擬采用軟件中的MODFLOW2005模塊進行三維流模擬。GMS地下水模擬系統(tǒng)是一個綜合性的、用于地下水模擬的圖形界面軟件。

1) 模型剖分

建模范圍為蘭考縣地熱開采區(qū)，總面積210 km2，垂向深度為2 000 m。模型水平方向剖分網(wǎng)格為22 514個，網(wǎng)格水平間距在150 m左右，開采井及回灌井位置進行網(wǎng)格局部加密以控制井損。綜合考慮地熱特征及地質條件，本次研究的3大熱儲層在垂向上剖分為5層：第1層為第四系含水層，底板埋深600 m，等效厚度200 m，忽略淺層水的影響；第2層為隔水層，底板埋深900 m，層厚300 m；第3層為明化鎮(zhèn)組熱儲層，底板埋深1 300 m，層厚400 m；第4層為隔水層，底板埋深1 800 m，層厚500 m；第5層為館陶組熱儲層，底板埋深2 000 m，層厚200 m。

2) 邊界條件

邊界條件是研究邊界處水位或流量應滿足的條件，反映模擬區(qū)內水流與周圍環(huán)境的相互制約關系。一般分為已知邊界水位變化規(guī)律的第一類邊界，已知邊界處流量變化規(guī)律的第二類邊界，已知邊界上水位與流量線性組合規(guī)律的第三類邊界。本次根據(jù)區(qū)內所有地熱井的實測水位數(shù)據(jù)，采用插值外推計算出邊界處的水位值。在可持續(xù)開采狀態(tài)下，應保持邊界處水位穩(wěn)定，即區(qū)內深層地熱開發(fā)不能對周邊區(qū)域熱儲層流場造成明顯影響，故將邊界定義為第一類邊界。

3) 初始水位

模型校正采用供暖末時刻同點位水位對比和供暖過程中動態(tài)監(jiān)測點水位對比相結合的方法。

4) 水文地質參數(shù)

考慮到第四系(Q)地下水經(jīng)過多年開發(fā)利用，其水位值下降劇烈，彈性釋水量已遠小于重力釋水量，故本次模擬中忽略其彈性釋水量。

3.2 初始賦值與模型修正

根據(jù)地熱井抽水試驗結果及監(jiān)測井長期實測數(shù)據(jù)對模型進行初始賦值。數(shù)值模型設定2018年底的地熱水流場為初始流場，第四系(Q)地下水初始流場如圖8a所示，區(qū)域內棕紅色越深表示水位埋深越大；明化鎮(zhèn)組熱儲層(N2)初始流場如圖8b所示，區(qū)域內藍色越深表示水位埋深越??；館陶組(N1)熱儲層初始流場如圖8c所示，區(qū)域內藍變綠越深表示水位埋深逐漸增大。

考慮到Q層和N2層開采井已經(jīng)被逐步關停，開采呈現(xiàn)出顯著減少的收斂趨勢，本次模擬Q層和N2層開采井只維持目前開采狀態(tài)(圖8a、圖8b)，重點關注層位為N1層，共布設長觀點5個，剔除異常數(shù)據(jù)后得到3個點位的有效動態(tài)監(jiān)測數(shù)據(jù)(圖8c)，3個點位分別為1、2、3號觀測點。為保證模擬符合實際，采用流場隨機識別法，即把各觀測孔的水位投放到相同時刻的模擬流場中的對應位置上，檢驗模擬效果并修正滲透系數(shù)等參數(shù)。通過模型擬合修正得到模型的運行參數(shù)值，詳見表3。水位擬合曲線如圖9—圖10所示。其中，圖9為空間上的擬合情況，圖10為時間上的擬合情況。實測值為各監(jiān)測井實際測得的水位值，修正值為通過參數(shù)調整后的模擬水位值。修正后模擬數(shù)值與實測值接近，變化趨勢保持一致，基本認定模型校正結果可以接受。

3.3 數(shù)值模擬結果

根據(jù)區(qū)內第四系(Q)、明化鎮(zhèn)組(N2)、館陶組(N1)3個熱儲層特征及地熱開采對地熱水初始流場的影響模式，綜合考慮時間離散等因素校正后的數(shù)值模型，每5天為一個步距，推演出不同時期熱儲層地熱水流場情況及其發(fā)展趨勢：采用限制條件下小規(guī)模開發(fā)利用方式模擬第四系(Q)及明化鎮(zhèn)組(N2)熱儲層在開采5 a(2023年)后的地熱水流場(圖11a—圖11b)；采用取熱后加壓同層回灌、封閉循環(huán)利用方式模擬館陶組(N1)在地熱開采5 a后的熱儲層熱水流場情況如圖11c所示。

圖8 不同熱儲層地下水初始流場Fig.8 Initial groundwater flow fields in different thermal reservoir

圖9 統(tǒng)測點水位擬合情況Fig.9 Water level fitting at survey points

圖10 1—3號動態(tài)監(jiān)測點水位擬合情況Fig.10 Water level fitting at monitoring point No.1-No.3

表3 水文地質參數(shù)Table 3 Hydrogeological parameters

3.4 數(shù)值模擬結果分析

通過對比蘭考縣深層地熱能開采至2023年后熱儲層流場的數(shù)值模擬結果與初始流場情況，可以看出：Q熱儲層流場由2處大面積降落漏斗狀態(tài)演變?yōu)樽晕飨驏|的面狀徑流，雖然Ⅲ區(qū)、Ⅳ區(qū)的水位顯著下降，但降落漏斗逐步消失；而N2地熱水流場漏斗面積小于0.05 km2，且未形成二級漏斗，流場漏斗效應明顯減小[19]；N1熱儲層流場漏斗效應明顯。

分析研究可知：Q熱儲層流場變化是由于其接受側向徑流補給，隨著該層位開采井的逐步關停，恢復到自西向東的原始流場模式；N2熱儲層隨著開采井分批關停，漏斗效應逐步減弱；N1熱儲層是區(qū)內地熱主要的開采層，隨著地熱資源持續(xù)開采，漏斗效應明顯，受開采井分布、熱儲層滲透系數(shù)、彈性釋水系數(shù)等影響，降落漏斗分布不規(guī)則。

4 結論

a.新近系館陶組(N1)地熱資源量明顯大于新近系明化鎮(zhèn)組(N2)資源量，回灌條件下前者地熱流體可采量為90 067 134.31 m3/a，后者地熱流體可開采量為16 468 660.80萬m3/a。

b.新近系明化鎮(zhèn)組(N2)熱儲和新近系館陶組(N1)熱儲的地質、地溫、水文地質、水化學特征明顯不同，為2個相互獨立的熱儲層。明化鎮(zhèn)組(N2)熱儲層的滲透系數(shù)與彈性釋水系數(shù)分別為1.00 m/d、0.001 73，館陶組(N1)的滲透系數(shù)與彈性釋水系數(shù)分別為1.25 m/d、0.047 20。

c.通過數(shù)值模擬，Q熱儲層流場隨著該層位開采井的逐步關停，恢復到自西向東的原始流場模式；N2熱儲層隨著開采井分批關停，漏斗效應逐步持續(xù)減弱；N1熱儲層隨著地熱資源持續(xù)開采，漏斗效應明顯，受開采井分布、熱儲層滲透系數(shù)、彈性釋水系數(shù)等影響，降落漏斗分布不規(guī)則。

d.采用加壓回灌方式可以大幅度提供回灌率，封閉循環(huán)利用地熱資源、取熱后同層回灌的方式能有效降低熱儲層流場的漏斗效應。具備條件時可同時開采古近系東營組的熱儲，綜合利用多層熱儲，合理規(guī)劃各層的開采量及平面布局，以有效控制漏斗效應。

圖11 2023年不同熱儲層地下水流場Fig.11 Groundwater flow field in different thermal reservoirs in 2023