干熱巖電學數值模擬方法與微觀特性分析

謝關寶

1 頁巖油氣富集機理與有效開發國家重點實驗室， 北京 102206 2 中石化石油工程技術研究院有限公司， 北京 102206

0 引言

干熱巖能源是一種高效低碳清潔能源，目前全球發現的干熱巖儲量是所有油氣總儲量的數十倍(冉恒謙等，2010；楊方等，2012；陸川和王貴玲，2015；許天福等，2012，2016；付亞榮等，2018)，其內部不存在流體或僅有少量流體(或致密不透水)的高溫巖體，儲層溫度一般在180 ℃以上，裂縫或斷層較為發育.作為一種重要的非常規資源(或能源)，研究其高溫條件下巖石物理性性質，可指導干熱巖熱儲資源的高效開發；但在一些特殊情況下，如裂縫巖石、高溫超高溫條件，常規巖石物理實驗難以開展，巖石物理響應及特征研究面臨極大挑戰.

近年來，數字巖心技術得到迅猛發展，通過構建數字巖心模型，分析巖石物理性質，已成為常規巖石物理實驗的有益補充.考慮巖石粒度分布，構建多孔隙數字巖心模型，可以精細表征巖石內部微觀結構，模擬實際巖石的導電特性(Bakke and ?ren，1997； ?ren and Bakke，2003).基于X-CT掃描構建碳酸鹽巖模型，可分析碳酸鹽巖力學特性(Arns et al.，2002，2004，2005a,b).對于頁巖及其他致密巖性，可以采用Nano-CT(Sisk et al.，2010；Song et al.，2018)、FIB-SEM(Sok et al.，2009；Curtis，2010；聶昕等，2016；Salvati and Korsunsky，2017)、MAPS(李潮流等，2016，王民等，2018)等多種技術手段，構建納米數字巖心，分析頁巖孔隙結構特征.干熱巖熱儲層位的巖石物理性質，在高溫高壓環境下會發生變化，尤其是高溫引起的電學性質變化，已成為識別干熱巖資源的重要特征之一(Spichak et al.，2007；Chen，2015).許多學者研究了電阻率與溫度之間的關系(Winsauer et al.，1952；Ucok et al.，1979；Kariya and Shankland，1983；孟銀生等，2010；余琪祥等，2011；童小龍等，2020).蘇聯科學家(Cheremensky，1982)，總結了巖石電阻率與溫度之間的經驗公式，但忽略了巖石熱開裂引起的電阻率異常變化.基于數字巖心分析技術，改變巖石微觀參數，可分析巖石的微觀物理性質(Adler et al.，1990；Blunt et al.，2002；Yue et al.，2011；Faisal et al.，2017；Tan et al.，2021).基于楓丹白露砂巖孔隙網絡模型，不考慮水濕巖石水膜傳導作用，在低含水飽和度情況下，電阻率指數的模擬結果遠大于實驗結果(Knackstedt et al.，2007)；如考慮水膜傳導作用，巖石電阻率模擬結果與實驗結果準確吻合(Liu et al.，2009).基于格子玻爾茲曼法，可模擬氣水分離過程，分析天然氣在地層水中的潤濕性和溶解度，分析其對儲層巖石電性的影響(Jiang et al.，2011).基于有限元法，采用消息傳遞機制(MPI)并行算法，可模擬超過2000×2000×2000個像素的巖石電性參數(Amabeoku et al.，2013).目前，基于數字巖心的巖石物理特性數值模擬研究面臨一個重大挑戰，即圖像分辨率與巖心尺寸(即圖像分辨率與儲層非均質性)之間的矛盾(孔強夫等，2015).

我國青海省共和盆地干熱巖發育，地熱資源豐富.本研究以該地區熱儲層巖心為研究對象，采用X-CT、QEMSCAN、Maps等多種實驗方法，構建了多尺度、多組分、高分辨率數字巖心；采用分形布朗運動方法(fBm)，構建了裂縫性數字巖心；考慮干熱巖儲層高溫的特點，在數值模擬中引入了溫度對巖石骨架的影響，并利用有限元法計算了巖石電阻率，分析了干熱巖電學特性微觀響應機理，為干熱巖熱儲層位評價和地熱資源開發提供了微觀理論基礎.

1 多尺度數字巖心構建

研究所用巖心取自共和盆地GR1井，其熱儲層位主要巖性為黑云母片麻巖、花崗巖和花崗閃長巖，孔隙度、滲透率都極低，裂縫和斷層較為發育(陳梓慧等，2015).研究鉆取了三顆巖心，如圖1所示，分別編號為No.1、No.2、No.3，對應取心深度為2980 m、3380 m和3600 m，巖性均為花崗巖.采用X射線衍射(XRD)及覆壓孔滲獲實驗取了巖心主要礦物成分與孔隙度等參數，如表1所示.實驗結果表明，研究區干熱巖主要造巖礦物為石英、斜長石和鉀長石，含少量黏土礦物；孔隙度非常小(最小為0.44%)，巖心較為致密，其中φ3>φ2>φ1.為了構建高分辨率數字巖心，采用QEMSCAN進行定量礦物分析，采用Maps進行微觀孔隙分析，采用X-CT掃描構建三維數字模型.

首先，對樣品進行二維大面積掃描電鏡成像(Maps，分辨率為250 nm，圖2a)和QEMSCAN礦物成分分析(圖2b)，確定巖心孔隙尺寸分布特征及礦物分布情況.根據巖心孔徑分布特征，選擇合適的分辨率進行X-CT掃描. 然后，根據礦物分布特征來確定納米級別的掃描區域(圖2b中的黃圈)，鉆取子樣品進行納米CT掃描，1號樣品子樣品記為GRY1，2號樣品子樣品記為GRY2，3號樣品子樣品記為GRY3-1.最后，由Maps圖像提取干熱巖樣品孔隙半徑特征，發現其孔隙尺寸很小，主要分布范圍為0.25～4.0 μm，后續采用的X-CT掃描分辨率為500 nm，其三維重構圖如圖2c所示，但此分辨率下孔隙識別度很低，孔隙度僅0.12%.因此，將X-CT掃描分辨率提高至65 nm，子樣品名稱為GRY3-2，其三維重構圖如圖2d所示，此時孔隙較為容易識別，孔隙度為1.39%.

圖1 干熱巖柱塞樣(a) 1號樣(25 mm); (b) 2號樣(25 mm); (c) 3號樣(25 mm).Fig.1 HDR plug cores(a) No.1 25 mm plug sample; (b) No.2 25 mm plug sample; (c) No.3 25 mm plug sample.

表1 干熱巖XRD和孔隙度實驗分析Table 1 XRD and porosity measurement experiment of HDR

Maps掃描區域與QEMSCAN掃描區域一致，將兩者圖像重疊對比，可以確定Maps圖像中每個像素點對應的灰度、礦物類型和微孔隙發育程度.根據Maps中礦物灰度值大致確定礦物密度排序，顏色越黑代表礦物密度越低，顏色越亮白代表礦物密度越高.通過3個Maps圖像分析可得到礦物的密度，密度大小排序為：ρ石英≤ρ斜長石<ρ鉀長石<ρ綠泥石+方解石<ρ黃鐵礦.這個密度排序為后續X-CT重建三維灰度圖像的多閾值分割提供了參考依據.

數字巖心組分劃分通常采用單閾值進行分割，0表示孔隙格架，1表示骨架，但單一閾值分割方式會喪失某些礦物成分信息.因此，在X-CT掃描圖像的基礎上，結合QEMSCAN礦物識別技術和XRD礦物組分分析，依據前述Maps獲取的礦物密度排序，對干熱巖數字巖心進行多閾值分割.數字巖心礦物組分分割結果顯示如圖2e、圖2f所示.可以看出，干熱巖樣品礦物以石英、斜長石和鉀長石為主，含少量綠泥石、蒙脫石等黏土礦物，含極少量鐵礦物，且巖心孔隙度均很小.從表2可以看到，從數字巖心得到的各礦物組分與XRD礦物分析結果接近，因此，可將數字巖心各礦物組分的電學綜合響應定義為骨架電阻率，根據理論方法和巖石物理實驗，構建骨架電阻率與溫度關系.通過多次掃描實驗構建的數字核心分別命名為GRS1、GRS2、GRS3-1、GRS3-2(65 nm).

2 電學特性微觀數值模擬方法

電學特性的數值模擬通常采用有限元方法(FEM)進行.考慮溫度對巖心電阻率的影響，將溫度對巖心電學特性的影響，轉換為其對巖心骨架的影響；將不同溫度下的骨架電導率作為有限元輸入，模擬不同溫度下整體數字巖心電流，得到整個數字巖心電阻率.

表2 干熱巖礦物成分對比分析Table 2 Companion of the mineral composition of HDR from XRD and constructed digital cores

圖2 多種掃描實驗及多組分高分辨三維數字巖心(a) Maps掃描結果； (b) QEMSCAN掃描結果； (c) 500 nmX-CT掃描結果； (d) 65 nmX-CT掃描結果； (e) 三維多組分模型； (f) 三維孔隙模型.Fig.2 Multiple scanning experiments and multi-component high-resolution 3D digital cores(a) High-resolution scanning by Maps; (b) mineral analysis by QEMSCAN; (c) Reconstruction of 3D grayscale model with 500 nm X-CT; (d) reconstruction of 3D grayscale model with 65 nm X-CT; (e) 3D multi component digital core, including potash (light purple), mixture of quartz and plagioclase (white), chlorite (red), and ferriferous minerals (yellow) ; (f) 3D pore space.

2.1 電阻率與溫度關系

干熱巖大多屬于酸性侵入巖，骨架導電能力很差，可看作絕緣體.當溫度T升高時，巖心電阻率ρcore會降低，巖石電阻率和溫度有如下關系(Cheremensky，1982)：

(1)

式中，α是經驗系數；ρ0初始電阻率(Ωm)；T0為初始溫度(℃)；ρcore是巖心電阻率(Ωm)；T為溫度(℃).

選取1號、2號、3號巖心子樣本測量其不同溫度下電阻率，樣品分別命名為No.1-1、No.1-2、No.1-3、No.2、No.3-1、No.3-2；測量樣品單軸加壓(12 MPa)，利用電阻絲對樣品加熱，在測量溫度點穩定15 min；測量溫度從室溫開始，至500 ℃，每50 ℃測量一次，每個測量點測量5次，取平均值作為該測量點數值，結果如圖3所示.電阻率的對數與溫度大致成線性相關特征，與式(1)的描述一致.擬合可以得式(1)中系數α(本地區為0.02)，ρ0在150 ℃時為7.0×108Ωm.

圖3 不同溫度干熱巖電阻率Fig.3 The measured HDR resistivity results at different temperatures

干熱巖體可視作一種多孔介質復合材料，包括骨架和孔隙兩部分.假定孔隙不導電，則巖心電阻率ρrock與骨架電阻率ρskeleton有如下關系式：

(2)

式中，φ為孔隙度(小數)；ρskeleton為骨架電阻率(Ωm)；ρrock為巖石電阻率(Ωm).

通過巖心測量的巖心電阻率，可求解骨架電阻率，將其代入有限元算法中，即可模擬不同溫度下干熱巖巖心電學參數.

2.2 有限元算法

有限元基本思想是變分原理，即將求解每個像素上的電壓分布問題，轉化為求解系統線性能量極值問題，從而確定數字巖心有效電導率.為使能量En取極小值，需滿足En對變量um(結點電壓)的偏導數均為零，即：

(3)

在三維數字巖心中，由于其自身就是由離散的像素構成，選定每一個像素為一個單元，單元的形狀為立方體，每個像素通過8個頂點與相鄰單元連接，像素的能量由各個節點上的電壓確定，則像素內的能量為

(4)

式中，σpq電導率張量；e(ex，ey，ez)是像素內點(x，y，z)處的局部場.

在有限元方法中每個像素的能量可以表示為

(5)

式中，矩陣Drs是三維模型中的8×8矩陣，也稱為“剛度矩陣”，ur是像素上第r(1，2,…，8)個節點處的電壓；us是像素上第s(1，2,…，8) 個節點處的電壓.

將所有像素能量相加，可以計算整個三維數字巖心的總能量，將系統總能量取最小值，計算所有節點電壓.如果已知每個像素單元的電導率張量，可以求得每個單元內的電流.然后，對三維數字巖心中所有單元電流取體積平均，得到三維數字巖心的等效電流.最后，根據電流與電壓的關系，得到巖樣等效電阻率.

3 電學微觀響應特征分析

設置模擬溫度10～500 ℃，每50 ℃模擬一次，數字巖心代表體積元設置為400體素，模擬結果以Z-Axis為準，計算不同溫度下的數字巖心電阻率，模擬結果如圖4所示.隨溫度升高，干熱巖電阻率呈減小變化趨勢，與式(1)的描述是一致的，且GRS1、GRS2數值模擬結果與實驗結果也相吻合；但GRS3-1、GRS3-2的數值模擬結果與實驗結果不一致，在其物理實驗中間段，變化范圍較大且呈現非線性特征，與式(1)相差較大.

是什么原因導致3號巖心電阻率變化趨勢與1號、2號不同呢？與1號、2號樣品相比，3號樣品埋藏深度更深，儲層溫度更高，熱應力超過巖心抗壓強度臨界值，導致熱開裂或巖心破裂，如圖5所示.此外，在3號樣品的Maps和X-CT圖像中可以清楚地看到一些微裂縫.實驗室加熱和原生地下高溫條件的共同作用導致3號樣品的裂縫更加發育，甚至破裂.

圖4 不同溫度下干熱巖電阻率數值模擬結果(a) GRS1； (b) GRS2； (c) GRS3-1； (d) GRS3-2.Fig.4 Resistivity numerical simulation results of HDR at different temperatures(a) GRS1； (b) GRS2； (c) GRS3-1； (d) GRS3-2.

圖5 3號巖樣裂縫顯示(a) Maps圖像； (b) 65 nm X-CT圖像； (c) 巖心開裂照片.Fig.5 Fractures in maps image; X-CT image and sample of No. 3 HDR(a) Maps image; (b) 65 nm X-CT image; (c) Rock sample after high temperature heating.

圖6為3號樣品的2個子樣品不同溫度下的電阻率數值模擬結果.由圖可見，隨著溫度升高，其電阻率在300 ℃后發生較大差異變化.圖6中虛線指示的是兩個子樣品的平均電阻率，當溫度達到300 ℃時，電阻率變化趨勢發生了變化，故推測此時巖樣產生了裂縫.因此，加熱過程中產生裂縫的溫度邊界約為300 ℃.

圖6 3號樣品不同溫度下電阻率： 熱開裂溫度邊界約為300 ℃Fig.6 Resistivity measurement results of two subsamples from the No.3 sample at different temperatures, showing the temperature boundary of thermal cracking is about 300 ℃.

為進一步分析熱裂解產生裂縫對電阻率的影響，研究采用分形布朗運動(fBm)方法在3號數字巖心中人為添加裂縫(Voss，1988； Kim， 2007)，形成帶裂縫數字巖心模型，命名為GRS3-3.因研究區裂縫多為低傾角，因此模擬中添加的裂縫均為水平裂縫.為了模擬加熱和熱開裂實際情況，模擬中隨著溫度的升高，生成的裂縫數量也相應增加.具體模擬過程及參數如下：當300 ℃時，添加第一條裂縫，其長度200體素，開度為2體素；當400 ℃時，添加第二條裂縫，450 ℃時添加第三條裂縫；隨著溫度升高，先前生成裂縫的尺度會增大(假設長度每升高50 ℃增長100體素)；當溫度升至450 ℃(或第一條裂縫貫穿巖心)時，則不再增加裂縫條數.整個添加裂縫的過程如圖7a—圖7c所示，最后構建的裂縫性數字巖心模型如圖7d所示.

對裂縫性數字巖心進行電學特性微觀數值模擬，結果如圖8所示.由圖可見，在300 ℃前，電阻率仍呈線性下降趨勢；加入裂縫后，由于裂縫內空氣不導電，電阻率減少趨勢變緩，并保持穩定(在巖心礦物相變前，如礦物脫水).此時，含裂縫數字巖心GRS3-3模擬的電阻率，與巖石實驗測得的平均電阻率變化趨勢一致，這說明裂縫性數字巖心構建與數值模擬方法是正確的.在獲取巖樣或進行高溫巖石物理實驗受限的情況下，利用上述數字巖心方法，能夠更方便地分析干熱巖高溫電學性質.

圖7 不同溫度增加裂縫模型示意圖(a) 350 ℃時裂縫模型； (b) 400 ℃時裂縫模型； (c) 450 ℃時裂縫模型 ； (d) 450 ℃時多裂縫模型.Fig.7 Schematic diagram of adding fractures at different temperatures(a) Adding the first fracture at 350 ℃；(b) Adding the second fracture at 400 ℃；(c) Adding the third fracture at 450 ℃； (d) The combined digital core with multipal fractures at 450 ℃.

圖8 3號巖樣電響應數值模擬結果Fig.8 Resistivity simulation results of two sub-cores from No.3 sample after adding fractures

4 結論

為分析不同溫度下干熱巖電學特性與響應機理，構建了多尺度、多組分、高分辨率數字巖心三維模型，分析了多尺度孔隙結構及和礦物成分特征，并研究了不同溫度下干熱巖電阻率變化規律.具體認識如下：

(1)采用X-CT、Maps、QEMSCAN和XRD多種實驗分析方法，構建了三維多尺度、多組分、高分辨率數字巖心模型，分析了干熱巖儲層孔隙特征和礦物成分特征，為研究干熱巖電學特性提供了有效載體.

(2)干熱巖電阻率隨溫度升高呈線性下降趨勢.對于深埋層伴有熱開裂的巖石，300 ℃前電阻率仍呈線性下降，300 ℃后電阻率變化幅度減小并趨于穩定.

(3)裂縫性數字巖心模擬與高溫巖石物理實驗結果一致，表明巖心加熱過程中可加熱過程中可能出現熱開裂.

(4)數字巖心技術可快速研究干熱巖高溫石物理特性，可有效解決現有巖石物理實驗無法定量分析的問題，為研究干熱巖熱儲層高溫巖石特性提供了一種新方法.

干熱巖高溫電學特性在300 ℃出現拐點現象，文中利用裂縫性數字巖心進行了驗證，物模結果與數模結果有較好的一致性；但溫度引起的干熱巖“裂紋到縫網再到破裂”這個過程是比較復雜的，還需綜合考慮干熱巖體的縫網拓展機理、巖心礦物相變、礦物會脫水、原位地層應力條件等因素，才能夠對干熱巖儲層巖石物理性質有更深刻、更準確的認識.