馬 峰,高 俊,王貴玲,3,劉桂宏,余鳴瀟,趙志宏,劉金俠

1.中國地質科學院水文地質環(huán)境地質研究所,石家莊 050061

2.自然資源部地熱與干熱巖勘查開發(fā)技術創(chuàng)新中心,石家莊 050061

3.中國地質科學院地球深部探測中心,北京 100037

4.清華大學土木工程系,北京 100084

5.中國石化集團新星石油有限責任公司,北京 100083

0 引言

地熱資源是21世紀最具開發(fā)潛力的可再生能源之一,可分為水熱型和干熱巖型[1-3]。雄安新區(qū)位于華北平原冀中坳陷中部,水熱型地熱資源豐富,其也是我國中東部地區(qū)地熱資源開發(fā)利用條件最好的地區(qū)之一[4],包括牛駝鎮(zhèn)地熱田、容城地熱田大部分區(qū)域以及高陽地熱田北部。雄安地熱以中、上元古界薊縣系碳酸鹽巖為主體熱儲,具有儲層厚度大、溫度高、滲透性好、水質優(yōu)、埋藏淺、易回灌的特點[5-6]。經(jīng)過多年的地熱開發(fā),雄縣地熱的利用已成規(guī)?；?是我國第一個用地熱替代燃煤供暖的“無煙城”[7],形成了政企合作、統(tǒng)一管理、集中回灌的“雄縣模式”[8]。其中,雄安新區(qū)牛駝鎮(zhèn)地熱田開發(fā)利用條件好,而容城地熱田相對于東部牛駝鎮(zhèn)地熱田而言,勘查程度較低。2017年以來,中國地質調查局加大了對容城地熱田的勘查力度,在容城凸起及周邊不同的構造位置部署地熱井10余眼,最大探測深度達到4 000 m,揭露了薊縣系高于莊組熱儲層,并對其產(chǎn)能進行了評價[9-10];2018年,在核心區(qū)鉆獲了容城地熱田產(chǎn)能最大的地熱井,單井供暖能力達到30萬m2[11];2019年,在啟動區(qū)鉆獲了水量112.8 m3/h、井口溫度71 ℃的高產(chǎn)能地熱井,并揭露了太古宙變質巖基底[4]。近年來,對容城地熱田的研究也逐漸增多,如:郭世炎等[12]對容城凸起地熱田儲層屬性與資源潛力進行了分析;馬峰等[4]對容城地熱田熱儲空間結構及資源潛力進行了分析,認為容城地熱田地熱資源賦存條件好,有利于地熱資源的規(guī)?；_發(fā)利用。薊縣系霧迷山組是雄安新區(qū)目前開發(fā)的主力熱儲層,其蓋層主要為新生界,其地層結構松散,導熱性差,具有良好的保溫隔熱作用[10]。

科學規(guī)范地開發(fā)利用地熱資源,可發(fā)揮地熱清潔能源的安全、穩(wěn)定、高效等特點。采灌均衡是將提取完熱量的尾水回灌到熱儲層,實現(xiàn)取熱不耗水,這是保持熱儲壓力、延長地熱田壽命、防止環(huán)境破壞及實現(xiàn)地熱資源可持續(xù)開發(fā)利用的有效模式[7,13-15]。地熱流體數(shù)值模擬普遍用于地熱流體資源的動態(tài)預測、開發(fā)設計、未來生產(chǎn)能力的評價及地下熱儲中熱水系統(tǒng)的運動規(guī)律研究等,已成為地熱流體開發(fā)和評價的重要技術方法[13-20]。大多數(shù)研究地熱回灌系統(tǒng)性質的解析方法都基于理想化的假設,比如含水層是水平的、均質的、各向同性的和等厚的等條件,在實際應用中具有很大的局限性。在現(xiàn)實的對井開采系統(tǒng)中,回灌和開采層可能位于熱儲層不同的深度,且具有不同的垂向和水平滲透性能。因此,采用數(shù)值方法解決這些問題具有很大的優(yōu)勢[21-22]。如:彭展翔[22]利用TOUGH2.0軟件建立天津霧迷山組熱儲的溫度場數(shù)值模型,研究不同采灌條件下熱儲層溫度場的變化趨勢,并對采灌方案進行了優(yōu)化;田光輝[23]利用GMS和TOUGH2.0軟件平臺分別對天津市溫泉度假村砂巖孔隙型熱儲和基巖巖溶裂隙型熱儲的可采量進行了評價;劉東林等[24]采用TOUGH2.0軟件平臺評價天津東麗湖地區(qū)霧迷山組地熱資源可開采潛力,計算得出了東麗湖地區(qū)霧迷山組可新增開采量;阮傳俠[25]建立了天津地區(qū)霧迷山組地質模型和三維數(shù)值模型,對其壓力場和溫度場進行了模擬計算和預測分析;陳建兵等[26]建立了地熱回灌井的溫度場和滲流場的三維數(shù)值模型,預測西安市三橋地區(qū)回灌井成井后井內(nèi)溫度場和滲流場的變化;宋美鈺等[27]利用TOUGH2.0 數(shù)值模擬軟件建立了地熱流體的數(shù)值模型,對天津地區(qū)霧迷山組熱儲壓力場變化趨勢進行了預測,分析得出地熱資源開發(fā)利用可采熱流體資源量與實際開采量的關系不大,而與回灌量有更為直接的聯(lián)系;王寬[28]采用TOUGH2.0 數(shù)值模擬軟件,根據(jù)研究區(qū)對井在預設情景下運轉30 a的模擬結果,對河南省延津縣新近系下部熱儲壓力場、溫度場進行模擬,數(shù)值模擬在淺層地熱能和干熱巖研究中得到了很好運用;馮波等[29]對U型地熱井供暖潛力進行了數(shù)值模擬,認為U型井式閉循環(huán)地熱系統(tǒng)水平井段長度為400～500 m,循環(huán)流速設置為80 m3/h,溫度為20 ℃左右時,可實現(xiàn)地熱能可持續(xù)開采;馬子涵等[30]建立了澳大利亞庫珀盆地三維分區(qū)均質滲透率模型和非均質滲透率模型,分別進行了熱儲溫度場、流場及采熱性能變化的研究。

目前雄安地區(qū)已開展了大量的數(shù)值模擬研究工作。例如:王樹芳等[31]開展了牛駝鎮(zhèn)地熱田回灌與示蹤試驗研究,建立了水平裂隙型介質模型,模擬結果顯示,400 m 以內(nèi)的井間距可能會導致開采井溫度快速下降,而 1 500 m 的井間距比較安全;劉久榮等[32]通過回灌試驗與數(shù)值模型,提出了雄縣牛駝鎮(zhèn)地熱田最優(yōu)開發(fā)利用模式為集中開采和集中回灌模式;龐菊梅[33]利用FEFLOW(finite element subsurface FLOW system)軟件對牛駝鎮(zhèn)地熱田南部地區(qū)進行水-熱耦合數(shù)值模擬,結果顯示,在現(xiàn)有的開采和回灌條件下,地熱田霧迷山組儲層開采井的地下水位呈下降趨勢,10 a內(nèi)地下水位下降約為 15 m,且在地熱水開采的中心區(qū)域出現(xiàn)地下水位降落漏斗;胡秋韻等[34]通過分析容城凸起地熱資源成藏模式,基于COMSOL軟件模擬了不同采灌量對可采地熱資源量的影響,通過在各開采區(qū)內(nèi)合理布置開采和回灌井,評價了容東地區(qū)采灌均衡條件下的地熱資源量。

雄安地區(qū)地熱資源的規(guī)?；?、可持續(xù)開發(fā)需要制定科學合理的采灌方案。本文基于COMSOL Multiphysics軟件,建立了容城地熱田地區(qū)整個薊縣系碳酸鹽巖熱儲地質模型和數(shù)值模型,對單井開采模式下薊縣系熱儲未來50 a內(nèi)溫度場和壓力場變化趨勢進行了預測;并針對下降趨勢制定了回灌方案,分析不同井間距對熱儲溫度場和壓力場的影響程度,以期為雄安地區(qū)地熱資源合理開發(fā)、可持續(xù)開發(fā)提供依據(jù)。

1 研究區(qū)地熱地質概況

雄安新區(qū)位于渤海灣盆地冀中坳陷北部,冀中坳陷內(nèi)部凹凸相間的構造格局和斷裂分布為熱量的傳導提供了有利條件[35-38]。研究區(qū)的凹陷區(qū)有巨厚的低熱導率沉積層,而凸起區(qū)碳酸鹽巖的高熱導率地層為熱流的側向運移提供了有利條件[39],從而導致了熱量向凸起區(qū)聚集,并在區(qū)內(nèi)形成了容城地熱田、牛駝鎮(zhèn)地熱田和高陽地熱田,本次研究的重點區(qū)域為容城地熱田(圖1)。

F1. 太行山山前斷裂;F2. 大興斷裂;F3. 容城斷裂;F4. 牛東斷裂;F5. 南苑--通州斷裂;F6. 大興斷裂;F7. 燕山山前斷裂。據(jù)文獻[38]修改。

容城地熱田主要地層包括新生界第四系、新近系,中元古界薊縣系、長城系和太古宇。薊縣系由鐵嶺組、洪水莊組、霧迷山組、楊莊組和高于莊組組成,其中霧迷山組為薊縣系目前開發(fā)的主力熱儲,在整個容城區(qū)域均有分布(圖2)。新生界是碳酸鹽巖熱儲的主要蓋層,其中第四系厚度為340～420 m,新近系明化鎮(zhèn)組厚度為400～500 m,新近系明化鎮(zhèn)組底界埋深在755～1 000 m之間。薊縣系厚度變化較大,為400～2 500 m,最大埋藏深度超過3 500 m。

圖2 雄安新區(qū)容城地熱田霧迷山組熱儲埋深等值線

2 容城地熱田數(shù)值模型

2.1 數(shù)值模擬方法

一般的數(shù)值模擬方法應用地熱井簡化模型。對于細長的地熱井結構,其長度與半徑之比可達 104以上,這是熱儲建模中一個極具挑戰(zhàn)性的數(shù)值問題。為了克服這一數(shù)值難題,Al-Khoury等[40-41]和Saeid等[42]提出了一種用一維線單元模擬地熱井中傳熱過程的偽三維井筒模型。本文進一步擴展了此模型,引入一維線單元對該三維井筒結構進行簡化,考慮沿井筒軸向的滲流傳熱過程,井筒內(nèi)流體與蓋層的熱交換通過等效換熱系數(shù)來近似考慮[34,43](圖3),其他主要假設條件[44]包括:1)儲層處于完全飽和狀態(tài),且不考慮儲層內(nèi)氣液相變過程;2)蓋層處于完全干燥狀態(tài),且不考慮不同儲層之間的水力聯(lián)系;3)井筒內(nèi)流體沿軸向流動,且同一深度流速沿徑向處處相等;4)考慮井筒內(nèi)流體性質如密度、黏度、熱傳導系數(shù)、比熱容等與溫度的相關性,但不考慮井筒內(nèi)氣液相變過程。

r1. 井徑;r0. 套管后井徑;r2. 灌漿后井徑;Δr. 井徑差;T0. 井壁溫度;T1. 套管井壁溫度;T2. 灌漿井壁溫度;Text.外部溫度;h. 深度。

本次建模采用成熟的地質體多場耦合有限元計算軟件COMSOL Multiphysics。COMSOL為一款通用的工程仿真軟件平臺,通過在一個模塊中添加其余模塊的因變量,實現(xiàn)各個模塊間的任意耦合。其優(yōu)勢在于用戶可以耦合多個物理場、化學場,并提供強大的開放接口供用戶修改控制方程。

2.2 模型范圍

數(shù)值模型計算范圍水平方向上為容城地熱田,面積約為300 km2;垂直方向上,根據(jù)巖性差異從上往下主要考慮第四系、新近系和薊縣系,根據(jù)容城已有地熱井的相關鉆井資料,地熱井井深在1 700 m左右,模型深度設為2 000 m。第四系和新近系的底板埋深根據(jù)已有的鉆孔資料,利用surfer軟件插值后導入COMSOL軟件中,建立的容城地熱田地質模型如圖4所示。

Q. 第四系;Nm. 明化鎮(zhèn)組;Jx. 薊縣系。

2.3 熱儲層參數(shù)

研究區(qū)薊縣系熱儲層參數(shù)按照表1取值。

表1 模型參數(shù)取值列表

深部熱儲中的流體密度(ρf)、黏度(μ)、導熱系數(shù)(kf)和比熱容(Cp,f)隨溫度(T)變化情況按如下經(jīng)驗公式[45-46]計算:

ρf=838.47+1.40T+0.003T2+3.72×10-7T3;

(1)

μ=1.38-0.02T+1.36×10-4T2-4.64×10-7T3+

8.90×10-10T4-9.08×10-13T5+3.85×10-16T6;

(2)

kf=-0.87+0.009T-1.58×10-5T2+7.98×10-9T3;

(3)

CP,f=12010.15+80.41T+0.31T2-5.38×10-4T3+

3.62×10-7T4。

(4)

2.4 網(wǎng)格剖分

在COMSOL Multiphysics中,綜合考慮計算精度與計算效率之間的平衡,地熱井及周圍區(qū)域的網(wǎng)格精細化,最小單元尺寸5 m,邊界及斷裂面網(wǎng)格細化,模型共包含約26萬個單元(圖5)。

圖5 雄安新區(qū)容城地熱田地質模型網(wǎng)格剖分示意圖

2.5 時間離散

由資料①可知,地熱井只在每年11月15日至翌年3月15日4個月供暖期內(nèi)處于開采或回灌運行狀態(tài),其余時間地熱井均處于停運狀態(tài),故地熱井流量可表示為

(5)

式中:q為地熱井流量,m3/h;n為開采周期數(shù);t為運行時間,月。

2.6 模型驗證

1)水位埋深擬合

容城區(qū)目前有兩眼地熱井的水位監(jiān)測數(shù)據(jù),分別是領秀城1井和領秀城2井。領秀城1井的水位監(jiān)測從2016-07-01—2017-10-30,領秀城2井的水位監(jiān)測則從2013-11-15—2017-10-30。通過調整模型的水文地質參數(shù)和相關邊界條件,使模擬水位與實測水位基本吻合(圖6),可利用擬合得到的參數(shù)進行容城區(qū)的采灌優(yōu)化設計。從水位分布云圖(圖7)可以看到,區(qū)內(nèi)因抽水而形成了漏斗,回灌則能有效地恢復儲層壓力。

a. 領秀城1井,2016-07-01—2017-10-30;b. 領秀城2井,2013-11-15—2017-10-30。

圖中黑色箭頭表示熱流運移方向。

2)溫度擬合

根據(jù)鉆井資料(1)河北雄安新區(qū)管理委員會綜合執(zhí)法局.雄安新區(qū)地熱資源保護與開發(fā)利用規(guī)劃(2019--2025年),雄安新區(qū):河北雄安新區(qū)管理委員會,2019.,利用surfer軟件進行插值后得到的模型地溫梯度如圖8所示,基巖地溫梯度基本在1.5 ℃/hm左右。領秀城1井的監(jiān)測水溫為51 ℃,領秀城2井的監(jiān)測水溫為49 ℃。通過調整相關的熱物性參數(shù),使得計算得到的兩眼監(jiān)測井的溫度與實測溫度一致,從而確定熱儲層參數(shù),如表2所示。

表2 研究區(qū)地層參數(shù)取值

圖8 研究區(qū)新生界地溫梯度

3 結果與討論

3.1 熱儲開采與回灌數(shù)值模擬

為控制霧迷山組熱儲層地熱水水位的持續(xù)下降和地熱尾水排放污染,地熱供暖尾水必須回灌,以保證采灌均衡。研究區(qū)采用“一采一灌”對井模式,目前已有11眼開采井和11眼回灌井,擬部署9對對井。按照圖9的井位分布設定,部署鉆孔的位置處也匹配1眼回灌井。井間距初步設計為600 m。開采量和回灌量均設置為100 m3/h,區(qū)內(nèi)已有地熱井的流量根據(jù)抽水試驗,按最大涌水量考慮,具體取值如表3所示。通過調整模型的水文地質參數(shù)和相關邊界進行模擬,結果如圖10、11所示。圖10為模擬不同開采年份的熱儲溫度分布。從圖10可以看出,隨著開采年限的增加,回灌冷鋒面在地下滲流場的作用下逐漸向開采井運移,從而導致某些開采井發(fā)生熱突破,開采溫度逐漸降低。本模型以開采井的平均溫度降低1 ℃作為熱突破的臨界值,不同開采井的溫度變化曲線如圖11所示。圖11中曲線的振蕩幅度與時間周期有關。統(tǒng)計開采井的熱突破時間可知:谷豐3井的熱突破時間最長,為25 a;谷豐2井的熱突破時間次之,為13 a,這是由于谷豐2井靠近回灌井,熱突破時間相對較短;領秀城1井的熱突破時間最短,約10 a左右,這是由于領秀城1井靠近回灌井領秀城4井,為避免快速發(fā)生熱突破可以采用減小采灌量的方法。

表3 研究區(qū)地熱井流量取值列表

圖9 研究區(qū)已有井位分布

a. 領秀城1井;b. 領秀城2井;c. 谷豐2井;d. 谷豐3井。

3.2 無回灌情況下

在沒有回灌井的情況下,對研究區(qū)已有的11眼開采井和擬部署的9眼開采井進行模擬。沒有回灌井的情況下開采井溫度通常不會降低,但長期的開采勢必會造成熱儲水位的下降,我們以熱儲層地熱流體開采最大允許埋深值150 m為約束條件,用來計算研究區(qū)地熱資源的最大開采量。模擬期從2017年開始,模擬周期50 a,計算結果如圖12所示,可以看到,在無回灌井開采40 a的情況下,在地熱井比較集中的區(qū)域水位埋深較大,并且因抽水形成的漏斗區(qū)域逐漸擴大,相對孤立的地熱井則會形成小范圍的抽水漏斗。通過監(jiān)測地熱井的水位埋深變化,可以得到地熱井的水位變化曲線(圖13)。從圖13可以看到,在40 a左右水位埋深降低至150 m,每年水位下降約1.13 m。

圖13 研究區(qū)開采井水位變化曲線

3.3 有回灌情況下

在有回灌的情況下,采用“一采一灌”方案的對井系統(tǒng),即1眼開采井匹配1眼回灌井,這樣熱儲水位由于回灌的作用,基本保持了壓力的平衡,唯一需要考慮的是在回灌條件下開采井的熱突破問題。我們以開采井的熱突破為限制條件,用來計算研究區(qū)地熱資源的最大開采量。計算過程為,在已有開采井和回灌井的基礎上(圖9),井間距分別考慮600、800和1 000 m 等3種情況,先按照圖14所示的布井方案在研究區(qū)域按最大數(shù)量布置對井系統(tǒng);然后計算對井系統(tǒng)開采井的熱突破時間,同時考慮在不同開采量下的開采溫度變化,計算時長為50 a。計算結果如圖15所示。最大開采量的統(tǒng)計結果如表4所示。

表4 研究區(qū)不同井間距對井系統(tǒng)的最大開采量計算

a. 井間距600 m;b. 井間距800 m;c. 井間距1 000 m。

a. 井間距600 m;b. 井間距800 m;c. 井間距1 000 m。

由表4可知:當井間距為600 m時,新區(qū)布井數(shù)量為138對,隨著開采量的增加(100、200到300 m3/h),熱突破時間從50 a降低到21 a,其中開采量為300 m3/h時,新區(qū)最大開采量為5.166 4×107m3/a,折合標準煤1.765 7×107t;當井間距為800 m時,新區(qū)布井數(shù)量為80對,隨著開采量的增加(300、400到500 m3/h),熱突破時間從51 a降低到28 a,其中開采量為500 m3/h時,新區(qū)最大開采量為5.257 4×107m3/a,折合標準煤1.796 8×107t;當井間距為1 000 m時,新區(qū)布井數(shù)量為50對,隨著開采量的增加(400、500到600 m3/h),熱突破時間從50 a降低到30 a,開采量為600 m3/h時,新區(qū)最大開采量為5.348 3×107m3/a,折合標準煤1.827 8×107t。計算結果表明:增大對井間距雖然減少了可布置的井對數(shù)量,但是延長了熱突破時間,增加了開采量,導致能開采更多的熱量。若以開采溫度降低1 ℃為熱突破標準,則在熱突破時間較長的情況下能開采更多的熱量。即在50 a 的開采年限里,開采溫度最終下降3和5 ℃ 的熱突破時間相差不多,最終開采的熱量相差不大(圖11)。

4 結論

1)通過COMSOL建立了容城地熱田群井回灌數(shù)值模型,根據(jù)已有的地質資料和相關的監(jiān)測數(shù)據(jù)對模型進行了識別和驗證,溫度和水位擬合結果較好,說明模型可以用于后續(xù)的計算。根據(jù)現(xiàn)有的采灌井布置方案,領秀城1井的熱突破時間最短,谷豐2井的熱突破時間較短,建議這兩眼井可以采用減小采灌量的方法,以延長熱突破時間。

2)采灌均衡情況下,采用“一采一灌”方案的對井系統(tǒng),考慮井間距分別為600、800和1 000 m 等3種情況,增大對井間距。雖然可布置的對井數(shù)量減少了,但是其開采量的增加及熱突破時間的變慢能開采更多的熱量,若以開采溫度降低1 ℃為熱突破標準,在開采量較小,熱突破時間較長的情況下能開采更多的熱量。在50 a的開采年限里,開采溫度最終下降3 ℃和5 ℃的熱突破時間相差不大,最終開采的熱量相差不多。建議在實際工程中,通過數(shù)值模擬方法選擇適宜的對井間距,保證開采量的增加及減慢熱突破時間,開采更多的熱量。