胡婷¹　芮振華^1,2

1.油氣資源與工程全國重點實驗室·中國石油大學（北京）

2.中國石油大學（北京）碳中和未來技術學院

0　引言

碳達峰、碳中和的時代背景下，二氧化碳捕集、利用與封存技術（CCUS）作為一種大規模的溫室氣體減排手段，成為助力“雙碳”目標實現的重要技術。地下廣袤的儲集空間提供了巨大的CO₂地質利用與封存潛力?？坍嫷叵聫碗s地質條件，并表征CO₂注入后的物理化學過程等相關技術水平的提高至關重要^[1]。注入油藏的CO₂一部分與原油發生相間傳質作用導致原油黏度降低，流動性增強，以提高原油采收率；另一部分CO₂溶解于地層水，甚至與巖石反應形成碳酸鹽礦物以實現安全封存。CO₂提高采收率與封存的本質是CO₂遷移與相態轉化的過程。因此，準確理解CO₂注入儲層后，CO₂的遷移轉化規律及分布相態演化，對CCUS工程的地質封存儲量評估以及安全性評價至關重要^[2]。

數值模擬技術是實現礦場碳封存潛力評價的主要方法，目前被廣泛應用的商業化碳封存模擬軟件是加拿大發布的CMG軟件[3]和美國發布的Eclipse軟件。此外，還有TOUGH、UNCONG、CO₂-SCREEN等其他模擬軟件。然而，由于沒有考慮地層化學反應的非均質性、未配套完善的地球化學反應參數數據庫而忽略可能形成的次生礦物，或只能實現用戶自行定義的簡單水—巖反應等原因，現有的成熟商業軟件不能構建符合實際地層條件的地球化學系統。實際地層中礦物種類復雜，除碳酸鹽礦物外，還有長石等硅酸鹽礦物和黏土礦物等鋁硅酸鹽礦物，而硅酸鹽、鋁硅酸鹽礦物的溶解正是為礦物態CO₂提供陽離子的主要途徑^[4]。礦物的溶解、沉淀反應大部分為反應動力學過程而非瞬間平衡過程，如不同酸堿度（pH值）的反應速率常數不同，影響礦物反應速率的反應表面積隨著礦物體積、潤濕性及離子半徑等參數變化而改變，礦物溶解及沉淀反應活化能參數的差異等，都構成礦物的復雜反應動力學過程^[5]。

目前的數值模擬軟件中，用戶自定義的簡單的礦化反應過程對于CO₂—水—油—巖系統復雜的地球化學反應過程的量化表征是不能完全實現的，對于地層的碳封存潛力，尤其是固相碳儲潛力的量化研究尚不完備。為此，通過程序開發手段，將閃蒸計算加入開源的反應溶質運移模擬軟件中，開發配套地球化學數據庫的熱—流—化耦合油藏CCUS數值模擬程序，利用改進后的軟件建立松遼盆地大情字井油田H59區塊的三維地質模型，厘清了CO₂—水—油—巖系統的相互作用機理，明確地球化學反應過程對滲流的影響機制，通過歷史擬合注采過程校正模型地層參數，量化表征不同注入階段及注采后CO₂遷移與相態轉化的時空演化過程，以期解決CO₂高效地質利用與封存的關鍵問題。

1 　基于熱—流—化耦合模擬軟件的開發

筆者基于熱—流—化耦合TOUGHREACT數值模擬軟件開展研究^[6]，該軟件是在原有的多相流體流動及熱傳導模擬軟件TOUGH^[7]的基礎框架上開發而成，軟件中的地球化學反應模塊與流體流動模塊為順序耦合^[8]，可考慮不同溫度、壓力、液體飽和度和離子強度條件下的地下熱物理化學過程^[9]。筆者在TOUGHREACT軟件的基礎上，增加了含多種烴類組分的油相，并將多種烴類組分的狀態方程計算和含油相的多相滲流計算耦合到該軟件的多相流體流動運算模塊中^[10]。

1.1 　相態數及主變量變化

對于儲層中可能存在的相態組合如表1所示。在該熱—流—化耦合數值模擬軟件的基礎上進行的模塊開發，盡可能地沿用了原有的主變量及其所在位置^[11]。若增加NHC個烴組分，加上水組分，那么整個系統的組分個數N_k=N_HC+1。N_HC個烴組分，需用（N_HC-1）個烴組分摩爾分數進行狀態刻畫。保留溫度（T）、壓力（p）參數及其位置，根據系統相態存在的情況主變量設置如表1所示。其表達式為：

式中z_i定義為某一相中烴組分i在總烴組分中的摩爾分數；m_HC,i 表示某一相中烴組分i的摩爾質量，mol；m_HC,k表示某一相中第k個烴組分的摩爾質量，mol；N_HC表示總烴組分數，個。

表1    相態組合及主變量設置表

保留原軟件中第一個相態位置（即序號為N_PH＝1）為氣相指針位置、第二個相態位置（即序號為N_PH＝2）為水相指針位置，設定油相的熱物性參數數組序號為N_PH＝3，代替原有NaCl析出形成的固相序號位置。因此，更改序列后的程序將原本N_PH＝3序列中的固相飽和度與孔隙度的變化相聯結。各相態滿足的關系式的變化如下式所示：

式中V_g、V_w、V_o分別表示地層中氣相、液相、油相態的體積，m³；?表示地層的孔隙度；V_m表示地層的體積，m³。

1.2 　狀態方程耦合多相滲流計算

針對系統中含有多組分烴的特征，利用PengRobinson（PR）立方方程^[12-13]對多組分在氣相或油相中的逸度進行計算，以此確定組分的相平衡狀態。該方程利用單組分方程結合二元相互作用系數計算多組分混合物，考慮了多組分間的相互作用，保證預測結果更準確，方程有2種表達形式，如式（3）、（4）所示。二元交互作用系數aij通常由實驗室測得的數據建立數據庫確定，其他如分子量、標態性質、偏心因子以及理想氣體的比熱系數等參數取自于已發布的屬性數據庫^[14]。

1.3 　耦合地球化學過程

在化學求解模塊中，所需求解的過程包括：①氣體組分CO₂溶解于水相形成HCO₃^–的過程；②水相中溶質遷移過程；③礦物溶解及沉淀過程。因此，在完成流動及熱流模塊的計算后，分別調用序列矩陣N_PH=1氣相中CO₂分壓及N_PH＝1氣相、N_PH＝2水相的流體流動速率及上一個時間步的地球化學系統的變化參數進行上述過程的計算。

如圖1所示，首先將考慮油相的流動及熱流計算模塊求解的流體流動速率導入溶質運移模塊，通過溶質運移計算確定各相態中的組分分布^[11]；之后，將氣相中CO₂分壓及水相中的各組分濃度分布導入化學計算模塊進行單網格的動力學平衡計算及收斂判斷；收斂平衡后，將更新的孔隙度（?）、絕對滲透率（K）及CO₂質量的變化量（）再傳回流動及熱流計算模塊，至此完成二次開發軟件中的流動及熱流計算模塊、化學結算模塊的順序耦合過程，實現對改進的流動及熱流計算模塊與原有的化學計算模塊的搭接^[18]。

圖1    CCUS-EOR數值模擬程序工作流程圖

1.4 　CO₂分布相態的計算方法

注入的CO₂最終變為4種不同的分布形態，即氣相、溶解于水相、溶解于油相、礦物固相。由于氣相的密度隨壓力變化大，因此，氣相CO₂封存量用質量守恒方法進行計算^[19]，各相態中CO₂分布量在改進的軟件中將通過以下公式進行計算：

2 　數值模擬及歷史擬合

2.1　模型構建

以中國CCUS工業化應用的示范場地——松遼盆地大情字井油田H59區塊為例，通過收集場地地層參數特征及井流監測數據，構建研究區的三維數值模型^[20]。目標層上白堊統嫩江組黑帝廟油層為輕質油，油密度為0.765 kg/m³。前人已開展的細管實驗表明，該油層的最小混相壓力為22.3 MPa，地層破裂壓力為49.0 MPa^[21-22]，因此，CO₂注入后可保證地層安全的條件下與原油發生混相^[23]。H59區塊采用正方形反九點井網進行開采，其目標儲層深度介于2 230～2 280 m，巖性以長石巖屑砂巖為主，平均孔隙度為12.7%，平均滲透率為3.0 mD^[24]。根據收集的井流數據資料，選取H59區塊邊緣位置的一個注采單元中的注入井H59-1井及開采井H59-3-3井建立注采模型，復現目標注采井13個月的歷史注采過程，并分析注采期內注入碳的遷移、轉化及分布相態特征，定量刻畫注采結束后第10年、第100年時儲層內注入CO₂的分布相態演化。

2.1.1　概念模型和邊界條件

目標注入井H59-1井及開采井H59-3-3井間距200 m，建立的單注單采三維場地數值模型不考慮H59-1井注采單元內其他開采井的影響。設置側向邊界位置距離注采井5 000 m及3 000 m，邊界網格設置體積無窮大，并設定為第一類邊界條件，降低邊界效應的影響。大情字井油田H59區塊黑帝廟油層單層厚度介于11.2～18.2 m^[25]，模型垂向厚度設置為12 m，井周網格相較于其他區域進行了加密剖分，加密部分網格數量為10×10×4共400個網格，每個網格在x、y、z方向上的距離分別為20 m、20 m、3 m（圖2）。概念模型設置在H59-1井儲層底部注入，在H59-3-3井儲層頂部開采，注入段和采出段儲層厚度均為3 m。

圖2    大情字井油田H59區塊概念模型圖

2.1.2　初始物性條件

根據大情字井H59區塊的地質資料，取各參數平均值為模型的初始物性條件，如表2所示。

表2    H59區塊黑帝廟油層初始物性條件數據表^[26]

初始原油擬組分依據前人研究結果設置，具體參數如表3所示。

表3    H59區塊黑帝廟油層初始原油擬組分數據表^[10]

相對滲透率參數^[19]設定如圖3所示，依據水—油及氣—油相對滲透率參數，利用STONE II模型計算氣—水—油三相系統中的流體相對滲透率。STONE II模型中相關參數設置如表4所示。

圖3    H59區塊黑帝廟油層儲層流體相對滲流率曲線圖

（資料來源：據本文參考文獻[19]）

表4    H59區塊黑帝廟油層STONE II模型參數設置表

2.1.3　初始化學條件

由于現場資料有限，模型中原始地層水的化學條件設定參考雷浩^[26]研究中的大情字井油田H59區塊的地層水化學組成進行設定，如表5所示。初始礦物組成設定為所收集資料中礦物含量的平均值^[27]，礦物組成及礦物的反應動力學參數如表6所示。

表5    H59區塊黑帝廟油層初始地層水離子組成及摩爾濃度數據表^[26]

表6    H59區塊黑帝廟油層初始礦物組成及其反應動力學參數數據表

2.2 H59區塊注采史擬合與預測

已開發軟件的可靠性和準確性已通過嚴格的對比驗證，并進行場地應用^[10]。針對目標區地質模型的參數設置進行歷史擬合，定量分析目標區CO₂的遷移與分布相態規律。

2009年5月至2010年5月，對注入井H59-1井井底壓力進行了連續監測，井流數據如圖4所示，圖中顯示了監測期內H59-1井注入水及CO₂的注入速率。H59-1井分3個階段進行注入，分別為注水階段、注CO₂階段和氣水同注階段，監測數據顯示井底壓力平均值為33.2 MPa。運用數值模型對監測的注入井井底壓力進行歷史擬合，并計算出開采井井底壓力為14.5 MPa。根據井流數據計算出向儲層中注入水和CO₂的累計量分別為0.654×10⁴ t、0.942×10⁴ t（圖5），依據井流數據復現H59-1井的流體注入過程，分析3種典型注入方式下的流體運移、轉化規律。

圖4    H59-1井生產曲線圖

（資料來源：據本文參考文獻[22]）

圖5    H59-1井的實際注入曲線圖

H59區塊H59-1井注采單元包括注入井H59-1井和開采井H59-3-3井，開采井H59-3-3井位于注入井H59-1井西南200 m處。數值模擬計算的生產數據如圖6所示，3個不同注入階段的日產油量、日產氣量均表現為增大的趨勢，在氣水同注時日產油量、日產氣量達到最大，而日產水量基本不變。這主要是由于注氣階段注入CO₂在注采井間大面積溶解于油相，改善了油相的流動性，有利于后續流體的注入及運移^[28]。此外，注采期3個階段共采出已注入的CO₂量為0.045×10⁴ t，主要來源于溶解在油相的CO₂，最終保留在儲層中的注入的CO₂量為0.897×10⁴ t。在注氣階段轉變為氣水同注階段時，模型計算得到的CO₂采出量迅速下降的現象與實測值變化情況一致，這說明軟件具有較好的場地應用效果。此外，在整體上計算值較實測值均勻分布，這主要是由于實際條件下生產數據受多種因素的影響，如地層非均質性、井間干擾等^[29-30]。

圖6    開采井H59-3-3井的模型計算產量曲線圖

（資料來源：實測數據來源于本文參考文獻[22]，有修改）

3 　注采期碳遷移和分布相態特征

3.1  CO₂遷移及空間分布特征

對注水階段、注氣階段及氣水同注階段地層含氣飽和度及含水飽和度的空間分布變化情況進行分析發現，在注采過程中，注氣階段及氣水同注階段注入的CO₂只在注入井端小范圍內呈現氣相，其余位置由于混相而以油相形式存在；含水飽和度的變化顯示，氣水同注階段注入水的運移范圍明顯大于注水階段，說明注CO₂能進一步增強后續注水的波及效率（圖7）。

圖7    數值模擬注采期注采井間流體空間分布剖面圖

注水前、注水階段、注氣階段及氣水同注階段含油飽和度的空間分布如圖8所示，儲層內含油飽和度在注水后顯著降低，由于在頂部開采，含油飽和度在不同深度上均呈現斜向上的蝶形分布。在注氣階段及氣水同注階段，由于CO₂與油相發生混相，靠近油相一側注入前緣的含油飽和度大面積增加。此外，溶解CO₂的油相密度降低而向儲層頂部擴散，驅動儲層頂部油相向下運移并與CO₂相互作用，不斷的傳質及擴散過程導致整個儲層厚度上的含油飽和度趨于一致^[31]。

圖8    數值模擬注采期注采井間含油飽和度剖面圖

注水階段、注氣階段及氣水同注階段油相中CO₂摩爾分數的空間分布如圖9所示，注氣階段氣相CO₂最遠遷移至距注入井160 m處。氣水同注階段CO₂已運移至開采井位置，由于氣、水的滲流速度存在差異，圖中局部低值顯示注入氣與注入水的滲流距離出現明顯差異現象，局部低值區域則是水段塞分布區。油相中CO₂摩爾分數以注入點為中心，呈現弧狀分布，這說明CO₂溶解油相后浮力作用減弱，油相中CO₂組分在垂向上趨于相對均勻分布^[32]。

圖9    數值模擬注采期注采井間油相中CO₂的摩爾分數剖面圖

3.2　注采期CO₂分布相態特征

確定地層中注入CO₂的分布相態是明確CO₂封存穩定性的關鍵，也是評估地層碳封存能力的前提。利用更新后的軟件對2009年5月至2010年5月間注入CO₂在氣相、溶解水相、溶解油相和礦物固相4種分布形態進行了定量表征（圖10）。由圖中可知，注采期內礦物以溶解為主，CO₂不存在以固相封存的相態。注氣階段CO₂主要溶解于油相，占比70.8%；其次為氣相，占比16.3%；溶解水相的CO₂占比12.9%。在氣水同注階段，氣相CO₂含量出現明顯拐點，該注入方式導致CO₂溶解水相的量明顯增加，同時氣相中CO₂含量降低?？紤]到溶解于水的CO₂能為礦物固相封存提供前提條件，因此，氣水同注是提高CO₂提采與安全封存雙贏的有效方式^[33]。氣水同注階段整體上占比最多的仍為溶解于油相的CO₂，占比72.2%；其次為溶解于水相的CO₂，占比17.8%；溶解氣相的CO₂占比10.0%。

圖10    數值模擬注入井H59-1井及開采井H59-3-3井注采期儲層內CO₂分布相態圖

4 　注采后碳遷移演化特征

4.1 CO₂分布相態演化

5　結論

1）研究區油藏CO₂混相驅條件下，CO₂在注入井端的小范圍內呈現氣相，在接觸到油相前緣后，CO₂會迅速與原油發生混相且受浮力作用減弱，油相中CO₂在垂向上趨于均勻分布，并向開采井端均勻推進。

2）氣水同注的方式導致CO₂溶解水相的含量明顯增加，為CO₂的穩定固相封存提供了條件，同時還保證了大部分注入的CO₂仍溶解于油相，這表明氣水同注的方式是實現CO₂提采與封存雙贏的有效方式。

3）注采結束時，注入的CO₂主要溶解于油相，占比72.2%；其次溶解于水相，占比17.8%；氣相占比10.0%，礦物固相可忽略不計；注采結束10年時，CO₂分布相態特征表現為溶解水相CO₂逐漸轉變為礦物固相，溶解在油相中的CO₂存在逐漸轉變為游離態氣相的趨勢，模型計算的單井注入量中礦物封存量占總量的0.01%，固相態封存量呈現穩定增長，其他相態含量均有所降低；注采結束100年時，地層中主要的沉淀礦物有方解石、鈉蒙脫石及菱鐵礦，主要的溶解礦物有鉀長石、鈣蒙脫石、鐵白云石及綠泥石。

4）大情字井油田H59區塊固碳礦物主要為方解石及菱鐵礦，主要礦化過程為鐵白云石及綠泥石溶解產生鐵離子、鈣離子與碳酸根結合形成固碳礦物方解石和菱鐵礦。