影響巖溶熱儲采出溫度的開采參數敏感性研究
——以河南清豐地區為例

韓 昀，李克文,*，何繼富，汪新偉，高楠安

(1.中國地質大學(北京) 能源學院，北京 100083；2.非常規天然氣地質評價與開發工程北京市重點實驗室，北京 100083；3.中石化新星(北京)新能源研究院有限公司，北京 100083；4.中國石化地熱資源開發利用重點實驗室，北京 100083)

地熱能是一種儲量豐富、分布廣泛、穩定可靠的可再生能源[1-2]。大力開發利用地熱能，對于實現“碳達峰碳中和”目標具有重要作用和深遠影響[3-4]。巖溶熱儲含水層被認為是火山地區之外最重要的熱水資源[5]，巖溶熱儲的采出溫度受到儲層特征參數的影響[6]。在構造裂縫或以溶洞為特征的地熱儲層背景下，裂縫長度和密度等裂縫屬性被認為是調節采出溫度的2個重要參數[7-8]。另外，斷層構造對采出溫度具有雙重影響。一方面，斷層可以作為熱水補給源的流動通道，提高采出溫度，增加地熱系統的可持續性和穩定性[9]。斷層通常會導致地溫梯度的變化，斷層處流體流速變化劇烈，從而改變地下溫度分布，并形成更高的溫度差，最終提高地熱能的開采速率和回收率。另一方面，斷層的存在可能引發地震和巖層位移的風險，增加地熱開發的地質風險和工程難度，還可能導致滲漏和流體損失的風險[10]。熱水可能通過裂隙和高滲透性區域流失，導致地熱能的浪費和采出溫度的降低。馬峰等[11]對河北雄安新區容城地熱田的碳酸鹽巖熱儲進行了研究，重點探究井間距對熱儲溫度場的影響。此外，其他熱儲回灌參數在很大程度上也決定著熱儲采出溫度。因此，對熱儲回灌參數的敏感性進行深入研究具有重要意義。

奧陶系巖溶熱儲是一種特殊的地熱儲層類型，具有高孔隙率、高滲透率、發育豐富的裂隙和溶洞等特征。以往研究已經對河南清豐奧陶系巖溶地熱儲層開發進行了一定的探索，主要圍繞地熱成因模式和地熱資源評價兩方面展開，如王迪等[12]對清豐地區巖溶熱儲進行了地熱地質特征的綜合調查。其研究成果顯示，清豐地區地熱田形成于正常大地熱流背景之下。這一地熱系統主要以太行山山脈和魯西隆起區的大氣降水為補給水源，經過地表滲透后沿著深大斷裂和巖溶不整合面滲透至地下，被認為是典型的中低溫傳導型地熱系統。趙婷[13]對清豐地區的熱源和地熱流體運移特征進行了重點分析，并明確了地熱的成因機理。一些學者[14-15]對渤海灣盆地內黃凸起地熱特征及其成因機制進行了系統研究。在地熱資源評價方面，章惠等[5]綜合鉆井數據與物探等資料，對渤海灣盆地內黃凸起東斜坡奧陶系巖溶熱儲進行了精細的地熱資源評價。但在明確斷層作用和回灌參數對采出溫度影響方面的研究仍然有限。鑒于此，筆者耦合達西定律、多孔介質傳熱和非等溫管道流物理場，建立考慮儲層非均質性條件下清豐地區的熱儲模型，并重點研究儲層地質條件、回灌速率、回灌溫度和注采井距等因素對采出溫度的影響，揭示區塊斷層對地熱開發的影響以及回灌參數影響巖溶熱儲采出溫度的作用機理，以期為巖溶地區地熱能開發奠定基礎。

1 清豐地區概況及幾何建模

1.1 研究區概況

研究區塊為地處渤海灣南側剪切構造區的清豐巖溶熱儲地熱藏，主要受黃驊?德州?東濮右旋橫推斷裂帶作用，具體位于內黃凸起背斜的南東側，西部和湯陰凹陷相接，北部為臨清坳陷，東部和東濮凹陷相鄰(圖1)。地質構造為單斜構造，地層向北東延伸，向南東稍有傾斜，地層傾角為3°～12°。內黃凸起是在4條大斷裂共同影響下形成[16]。

圖1 清豐地區構造位置Fig.1 Tectonic location of karst geothermal reservoirs in the Qingfeng area

地層按新老順序依次為奧陶系、石炭?二疊系、侏羅?白堊系、新近系和第四系，熱儲層以奧陶系臺地相碳酸鹽巖為主，地層分層及巖性特征見表1。

表1 地層分層及巖性Table 1 Strata and their lithologies

1.2 幾何建模

COMSOL Multiphysics軟件基于有限元方法，通過求解偏微分方程組來模擬真實的物理場景，其在多物理場耦合、建模以及非等溫流體流動模擬等方面擁有顯著優勢，是地熱開發仿真與模擬研究的重要工具[17]。首先根據清豐地區已有的7口地熱井測井數據(圖2)，利用Surfer軟件，采用克里金插值法繪制了清豐地區不同地層的地溫分布圖，并導入COMSOL軟件，結果如圖3所示。

圖2 清豐地區部分井深度與測井溫度關系Fig.2 Relationships between the depths and log-derived temperatures of some wells in the Qingfeng area

圖3 河南清豐地區不同地層頂面地溫分布Fig.3 Distribution of geotemperatures at top surfaces of different strata in the Qingfeng area,Henan Province

模型的孔隙率φ和滲透率k主要是利用Petrel軟件中的屬性建模模塊建立的。原始數據源自地熱井的測井數據，采用克里金插值法進行插值。處理后的數據分布如圖4所示。

圖4 插值后的滲透率和孔隙率分布直方圖Fig.4 Histograms showing the distributions of permeability and porosity after interpolation

其次，根據研究區熱儲的構造形態和儲層分布規律，選用COMSOL Multiphysics軟件，結合7口地熱井的井位坐標、斷層分布和地層分層數據，建立了儲層的幾何模型。利用Surfer軟件和地層分層數據，對各個地層界面埋藏深度進行克里金插值，形成各地層界面埋藏深度的插值數據。然后，將形成的數據導入COMSOL構建地層界面，并最終生成地層底面參數化曲面(圖5)。之后，根據本區巖溶熱儲中斷層和井位的分布情況，在模型中添加斷層面和地熱井(圖6)。通過COMSOL內置的幾何操作完成清豐地熱田的三維幾何建模(圖7)。模型大小為3 200 m×1 800 m×2 100 m。模型頂面設為定溫邊界，取年平均氣溫16℃，側面設為溫度開放邊界和水位邊界，底面設為隔水邊界和熱通量邊界，熱通量取當地大地熱流65 mW/m2[18-19]。

圖5 模型地層底面參數化曲面Fig.5 Parameterized bottoms of strata in the model

圖6 模型井位分布Fig.6 Well distribution in the model

圖7 三維地熱模型Fig.7 3D geothermal reservoir model

2 清豐地區巖溶熱儲數值模擬

2.1 模型假設和數學方程

2.1.1 模型假設

(1) 儲層巖石被視為由巖石基質和離散裂縫組成的三維裂縫多孔介質。巖石基質可簡化為連續多孔介質，巖石基質滲透率遠遠低于裂縫的滲透率，裂縫是儲層的主要流動通道；(2) 巖石基質和裂縫中的液體流動服從達西定律，不考慮流體流動的黏性效應；(3) 地熱采出井中流體流動速度斷面充分發展，速度不會在同一斷面處發生變化，井筒四周溫度相同；(4) 水和巖石基質之間的熱交換是通過對流和傳導過程實現的；(5)地層中的多孔介質熱對流和熱傳導控制的幾何域由多個固定成分的單個流體流經的多孔介質域構成[20]。

2.1.2 數學方程

多孔介質中流體流動連續性方程和采出井中非等溫管道流連續性方程[21]分別為：

多孔介質中流體流動的運動方程和采出井中非等溫管道流運動方程[22]分別為：

多孔介質中的傳熱控制方程[21]為：

多孔介質傳熱方程(5)中的第一項是滲流速度相關項，方程(3)為流體流動達西定律方程，其中速度與流體性質有關。通過聯立求解以上方程，得到模型的數學控制方程。

2.2 模型主要參數

考慮儲層非均質性，模型的地層參數見表2，模型的斷層參數見表3，孔隙率采用插值后的奧陶系儲層孔隙率。

表2 模型地層參數(巖石物理性質)Table 2 Parameters of strata in the model (rock physical properties)

表3 模型斷層參數Table 3 Fault parameters of the model

考慮到熱儲可持續開發的特征，采用同層對井回灌進行模擬，注入井的注水層與采出井的取水層均位于奧陶系，與實際現場注采井所處層位和井位一致。模擬時間為60 a，模型的生產參數見表4。

在COMSOL Multiphysics中使用自由四面體網格對整個幾何進行了剖分，剖分后的三維熱儲模型包含144 809個四面體單元、26 102個網格頂點和1 395個邊單元。其次，為了提高網格收斂度和模型模擬結果的準確性，將包含斷層的奧陶系地層網格進行了超細化處理。

2.3 模型驗證

為驗證建立的熱流耦合模型準確性，本文將模擬得到的采出溫度與清豐地區采出井實際現場監測的抽水溫度進行了對比。對比的時間范圍為實際現場抽水試驗時間，總計12 h，見表5。

表5 抽水溫度對比Table 5 Comparison of temperatures of pumped water

當抽水時間小于3 h時，實際生產井抽水溫度略高于模型生產井抽水溫度；當抽水時間大于3 h時，模型生產井抽水溫度開始高于實際生產井抽水溫度。為了更直觀地顯示模型的可靠性，將實際抽水與模型模擬抽水的對比溫度隨時間變化的情況繪制成柱狀圖(圖8)。從圖中可以看出，二者隨時間變化的趨勢基本一致，實際采出井抽水溫度略低于模擬溫度，平均相對誤差5.41%，誤差在合理范圍內，在一定程度上說明所建清豐地區的數值模擬模型是比較準確可靠的。

圖8 抽水溫度擬合Fig.8 Fitted temperatures of pumped water

3 斷層對熱儲采出溫度的影響

為了明確清豐地區斷層在熱儲采出溫度中的作用機制，從斷層影響儲層熱交換和熱儲采出溫度兩方面進行分析。

3.1 熱交換

圖9對比了開發60 a后有無斷層情況下的流體流速分布。從圖中可以看出，斷層的存在會改變流體流場分布，表現為斷層的存在導致了流體流速異常，整體表現出流速減小的趨勢。另外，在斷層附近，流體流速可能明顯增加或減小。導水性斷層會導致流體流速增加，阻水性斷層會導致流體流速減小。該區塊熱儲共發育4條斷層，根據模擬結果，4條斷層中2條為NNE向導水斷層，2條為NWW向阻水斷層。NNE向導水斷層可以增加地熱儲層滲透性，從而提高流體流動性和傳熱效率。斷層周圍的裂隙可以作為運移地熱流體的通道，增強熱儲巖石和流體之間的熱交換，而NWW向阻水斷層作用相反。從圖9a可以看出，盡管模型未考慮斷層的存在，但由于儲層非均質性的作用，無斷層作用下的地熱流體流速呈現出不均勻擴散形式。從圖9b可以看出，2條斷層中部的流體流速較其他地區要高。高速流體集中在采出井和注入井底部附近，流線方向相反。4條斷層圍成的區域地熱流體流速最高。

圖9 有無斷層情況下流速分布對比Fig.9 Comparison of flow-velocity distribution with or without faults

3.2 采出溫度

在開始取熱階段，較低溫度區域主要聚集在注入井附近。隨著開采時間的增加，這些較低溫度區域逐漸擴散至生產井，導致生產井中溫度急劇下降[25]，這種現象通常被稱為熱突破[26]。本文將采出溫度下降2℃的時間定義為熱突破時間。圖10為有無斷層情況下的熱儲采出溫度變化對比圖。從圖中可以看出，在開采前期，無斷層情況下采出溫度高于有斷層存在的情況。無斷層情況下的熱突破時間為36 a，有斷層情況下的熱突破時間為32 a，表明斷層的存在加速了熱儲中較低溫度區域向生產井的擴散，導致熱突破時間提前。

圖10 采出溫度變化對比Fig.10 Comparison of variations in temperature of produced geothermal water

在開采后期，當發生熱突破后，有斷層時的采出溫度遠高于無斷層情況，采出溫度隨時間的變化曲線斜率遠高于無斷層的情況，溫度下降幅度逐漸趨于平緩。表明斷層有助于提高清豐地區熱平衡時的采出溫度，最終有利于增加采熱量。原因可能在于隨著開采時間的推移，斷層可能會形成一些通道或通道網絡，為熱儲提供了熱量補給通道，通過斷層帶，高溫流體可以從深部地熱系統進入采出井，導致有斷層情況下的溫度下降變緩。

為了明確斷層提高熱儲采出溫度的作用機理，繪制了模擬60 a后有無斷層的溫度等值面分布對比圖(圖11)和熱通量運移對比圖(圖12)。根據圖11a和圖11c，采出井和注入井之間的NNE向導水斷層加速了注入流體沿斷層的擴散運移，因而延緩了冷鋒面向開采井的擴展。根據圖12a和圖12b，有無斷層的地層整體熱通量運移差別不大，原因是斷層對于熱通量的改變與井作用相比是微弱的，注入井不斷注入冷流體和開采井持續開采熱流體對于地層熱通量的改變起主導作用，但從圖12b中的藍色方框部分可以看出，存在由深部地熱系統通過斷層進入熱儲的熱通量，這表明斷層在一定程度上為熱儲提供了熱水補給源，通過斷層帶，高溫地熱流體可以從深部地熱系統進入采出井，增加地熱系統的可持續性和穩定性，表明該區塊斷層的存在增加了地熱能的利用潛力。

圖12 有無斷層情況下熱通量運移對比Fig.12 Comparison of heat flux transport with and without faults

4 熱儲回灌參數的敏感性分析

熱儲回灌參數在一定程度上控制著采出溫度，為了明確回灌參數對熱儲采出溫度的影響，分別對回灌溫度、回灌速率和注采井距進行了敏感性分析，并給出考慮儲層非均質性優化后的參數組合。

4.1 回灌溫度

圖13為模擬60 a后不同回灌溫度下儲層溫度分布。從圖中可以看出，不同回灌溫度下的儲層溫度呈現出一致的變化規律，均從注入井井底開始降低，并沿斷層向采出井突破。儲層東側的NNE向導水斷層的溫度變化幅度大于西側斷層，一方面，由于清豐地區奧陶系巖溶熱儲東側的NNE向斷層的導水性在一定程度上優于西側斷層，另一方面布井方式也會影響儲層溫度分布。

圖13 不同回灌溫度下儲層溫度分布Fig.13 Distribution of geothermal reservoir temperatures under different reinjection temperatures

圖14為模擬60 a后不同回灌溫度下采出溫度隨時間的變化。從圖中可以看出，不同回灌溫度下采出溫度變化不大，并呈現一致的變化規律，回灌溫度為10℃的情況下，熱儲采出溫度最低，隨著開采時間的增加，熱儲開發在高回灌溫度下的優勢逐漸顯現，溫度越高越有利于維持熱儲系統的穩定性和可持續性，并保持較高的地熱能利用效率。隨著熱儲開發進入中后期，熱儲熱能的賦存量逐漸減少，回灌溫度相對較高有助于減緩這一過程，使熱儲系統的運行時間更長。然而，在回灌溫度為20℃的情況下，熱儲系統可能已經實現與周圍環境的熱平衡，因此，再提高回灌溫度對熱儲采出溫度的影響不大，考慮到地熱持續開發的經濟成本，20℃為清豐地區奧陶系巖溶熱儲的最優回灌溫度，此發現為該地區熱儲開發方案的優化提供了重要參考。

圖14 不同回灌溫度下采出溫度隨時間變化Fig.14 Time-varying temperature of produced geothermal water under different reinjection temperatures

4.2 回灌速率

圖15為模擬60 a后不同回灌速率下儲層溫度分布。從圖中可以看出，較高的回灌速率會明顯加快熱儲溫度的下降并縮短熱突破的時間，這是因為大量冷流體會吸收儲層中的熱量，導致熱儲溫度降低。此外，高回灌速率加快了地熱流體的運移速率，降低了熱儲巖石和流體之間的熱交換時間。因此，過高的回灌速率不利于地熱能的可持續開發利用。

圖15 不同回灌速率下儲層溫度分布Fig.15 Distribution of geothermal reservoirs temperature under different reinjection rates

為得到準確的最優回灌速率，在不考慮熱儲水位下降的前提下，擴大回灌速率參數的研究范圍，繪制了模擬60 a后不同回灌速率下采出溫度隨時間的變化圖(圖16)。從圖中可以看出，開采40 a后，以50 kg/s的回灌速率為分界點，高于50 kg/s回灌速率下的采出溫度隨開采時間快速下降，開采60 a高回灌速率下曲線的末端斜率遠大于回灌速率低于50 kg/s的情況，這意味著其開采溫度將隨開采時間持續降低?；毓嗨俾实陀?0 kg/s時，隨著回灌速率的增大，采出溫度逐漸增大?；毓嗨俾矢哂?0 kg/s時，隨著回灌速率的增大，采出溫度逐漸降低。原因在于當回灌速率較低時，地熱水沒有得到足夠的有效補充，在此情況下，采出井采出流體為儲層中原本溫度較高的流體。高回灌速率雖然能夠及時補充開采井開采導致的熱儲流體損耗，保持儲層壓力，但由于回灌溫度一般低于采出溫度，大量冷水注入會導致采出溫度的降低。

圖16 不同回灌速率下采出溫度隨時間變化Fig.16 Time-varying temperature of produced geothermal water under different reinjection rates

綜合來看，50 kg/s為最優回灌速率。為保持地下熱水資源的可持續利用，合理調整回灌速率是必要的。

4.3 注采井距

圖17為模擬60 a后不同注采井距下儲層溫度分布。從圖中可以看出，當注采井距為505 m時，熱儲溫度變化最劇烈，儲層溫度降低幅度明顯，導致儲層更快地熱突破。隨著井距的增加，冷鋒面到達生產井的時間逐漸延長。當注采井距為1 000 m時，由于熱儲南部NNE向導水斷層的存在，加快了注入流體在斷層的擴散，導致斷層附近熱儲溫度明顯下降。另外，較大的注采井距會導致儲層中形成更大的溫度梯度，熱量需要通過較長的距離傳導才能從采出井采出。并且，大注采井距還會導致儲層溫度分布在儲層厚度方向上發生變化。

圖17 不同注采井距下儲層溫度分布Fig.17 Distribution of geothermal reservoir temperatures under different well spacings

圖18為模擬60 a后不同注采井距下采出溫度隨時間的變化。從圖中可以看出，采出溫度隨開采時間的變化規律以注采井距340 m為分界，當注采井距小于340 m時，采出溫度在開采初期劇烈下降，隨后趨于平緩；當注采井距大于340 m時，注采井距的增加呈現出一致的變化趨勢，當注采井距340 m時，采出溫度在開采60 a的降幅達到15.8℃，這是由于過小的井距導致了更快時間的熱突破。理論上，在滿足儲層壓力下降的前提下，較大的注采井距意味著流體在儲層中停留的時間更長，增加流體與儲層的換熱時間，從而降低熱儲由于回灌作用導致的溫度下降程度。但在此研究中，觀察到清豐地區的巖溶熱儲開發效果受到斷層和儲層非均質性的雙重影響。這些因素在不同注采井距下對采出溫度產生顯著影響。具體而言，同樣是滿足熱儲壓力下降條件，當注采井距835 m時，采出溫度高于注采井距1 000 m的情況。原因在于井距1 000 m時，注入井與儲層西側NNE向導水斷層距離過近，造成注入流體更快地擴散，在儲層非均質性作用下，儲層內部的溫度變化較大，采出溫度出現低于注采井距835 m的情況。

圖18 不同注采井距下采出溫度隨時間變化Fig.18 Time-varying temperature of produced geothermal water under different well spacings

該結果說明了地熱開發中斷層和儲層非均質性的重要性，并揭示不同注采井距下采出溫度的差異。這表明熱儲回灌參數的優化需要考慮儲層非均質性，忽略斷層的存在或將儲層視為均質，在一定程度上會導致參數優化的偏差。在地熱開發中，需要充分考慮這些因素，以制定優化的注采方案，實現更高效、可持續地利用地熱能。

綜合以上研究結果，給出了優化后的清豐地區特定場景下的熱儲開發主要參數：回灌溫度為20℃、回灌速率為50 kg/s、注采井距為835 m。

5 結論

a.在高于50 kg/s的回灌速率下，隨著開采時間的增加，采出溫度呈現快速下降的趨勢。盡管高回灌速率有助于及時補充熱儲水量的損耗，維持儲層壓力，但由于回灌溫度較低，大量冷水注入導致采出溫度降低。因此，在實際開采中應選擇合適的回灌速率，避免過高速率導致溫度下降過快的問題。

b.在開采前期，無斷層情況下采出溫度高于有斷層情況，這是由于斷層的存在加速了熱儲中低溫區向生產井的擴散，導致熱突破時間提前。然而，在開采后期，有斷層時的采出溫度遠高于無斷層情況。

c.在儲層非均質和斷層的影響下，注采井距835 m時的采出溫度高于1 000 m的情況，說明儲層非均質性和斷層對開采參數的優化有重要影響。

d.基于考慮儲層非均質性的數值模擬結果，優化了清豐地區特定場景下的回灌參數：回灌溫度為20℃、回灌速率為50 kg/s、注采井距為835 m。

e.未來應考慮儲層沉降等復雜條件下建立更加符合實際情況的熱儲模型，進一步提高數值模擬結果的準確度，有利于制定更加切實可行的熱儲開發方案。

符號注釋：

A為采出井井筒斷面面積，m2；Cp為常壓下流體比熱容，J/(kg·K)；dh為平均水力直徑，m；fD為達西摩擦因子；F為體積力，N；g為重力加速度，m/s2；k為多孔介質滲透率，m2；p為水頭壓力，Pa；q為熱通量，W/m2；Q為熱源，W/m3；Qm為流量，kg/(m3·s)；t為時間，s；u為流體流速，m/s；z為深度，m；?為哈密頓算子；ρ為流體密度，kg/m3；(ρCp)eff為多孔介質的等效體積熱容，J/(cm3·K)；φ為巖石孔隙率；μ為流體黏度，Pa·s；λl為流體導熱系數，W/(m·K)；λd為各向導熱系數，W/(m·K)；λe為有效導熱系數，W/(m·K)；λp為多孔基質導熱系數，W/(m·K)；θp為多孔基質體積分數，%。