雄安新區容城地熱田碳酸鹽巖熱儲采灌數值模擬

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

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

2.自然資源部地熱與干熱巖勘查開發技術創新中心,石家莊 050061

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

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

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

0 引言

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

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

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

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

1 研究區地熱地質概況

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

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

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

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

2 容城地熱田數值模型

2.1 數值模擬方法

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

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

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

2.2 模型范圍

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

Q. 第四系;Nm. 明化鎮組;Jx. 薊縣系。

2.3 熱儲層參數

研究區薊縣系熱儲層參數按照表1取值。

表1 模型參數取值列表

深部熱儲中的流體密度(ρf)、黏度(μ)、導熱系數(kf)和比熱容(Cp,f)隨溫度(T)變化情況按如下經驗公式[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 網格剖分

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

圖5 雄安新區容城地熱田地質模型網格剖分示意圖

2.5 時間離散

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

(5)

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

2.6 模型驗證

1)水位埋深擬合

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

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

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

2)溫度擬合

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

表2 研究區地層參數取值

圖8 研究區新生界地溫梯度

3 結果與討論

3.1 熱儲開采與回灌數值模擬

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

表3 研究區地熱井流量取值列表

圖9 研究區已有井位分布

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

3.2 無回灌情況下

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

圖13 研究區開采井水位變化曲線

3.3 有回灌情況下

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

表4 研究區不同井間距對井系統的最大開采量計算

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

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

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

4 結論

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

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