二氧化碳地質封存中蓋層力學完整性數值模擬研究綜述

劉苗苗，孟令東，王海學，吳 桐

(1.東北石油大學，河北 秦皇島 066004；2.黑龍江省油氣藏及地下儲庫完整性評價重點實驗室，黑龍江 大慶 163318；3.黑龍江省石油大數據與智能分析重點實驗室，黑龍江 大慶 163318)

0 引 言

CO2地質封存是控制大氣中CO2濃度最科學有效的手段之一[1]。將大量CO2持續注入地質封存體內，會引發一定范圍內的壓力傳播、流體運移等現象[2-3]，對蓋層力學穩定性產生重要影響[4]。要使CO2持續注入并穩定封存幾百或幾千年不逸散，就要保證地層壓力不超過蓋層突破壓力，因此，需定量表征封存過程中蓋層巖石的變形損傷程度，正確評價其力學完整性[5]。針對該問題，目前的研究方法主要有理論分析法[6-7]、室內實驗[8]或場地淺層含水層的注氣試驗法[9]、數值模擬法[10-14]。相比之下，數值模擬法能有效彌補其他研究手段的缺陷，已成為當前地質工程學中定量化研究的有效方法之一。鑒于此，對CO2地質封存中蓋層力學完整性分析相關問題的數值模擬研究現狀進行了綜合分析，以期為CO2地質封存的相關研究提供技術支持與引導。

1 CO2地質封存研究中的數值模擬方法

大量CO2的持續注入會導致巖石所處的應力狀態、流體的滲流路徑等發生變化[15]，對這一過程的刻畫是分析巖層應力-應變規律的關鍵。室內實驗方法只能分析該過程產生的機理，而數值模擬方法能充分考慮各個過程間復雜的耦合關系，可直觀有效地對地質模型進行數學和力學分析。目前，已有許多學者從不同角度，使用不同的軟件、模型及方法對CO2地質封存的相關問題進行了數值模擬研究。

1.1 多場耦合原理

研究顯示，將CO2以超臨界態注入咸水層封存是最為理想的方式，而CO2咸水層封存則是一個復雜的熱力學-流體力學-巖石力學-化學多場相互作用的過程[16]，如圖1所示。由于CO2注入過程一般只有幾十年，這期間化學場對蓋層力學完整性的影響作用幾乎可以忽略。此外，研究顯示，CO2注入引起的溫度變化一般僅局限于注入井周圍的幾十米[17]。因此，許多蓋層力學完整性數值模擬研究均忽略溫度的變化，僅研究滲流-應力耦合過程。根據前人研究[18-19]，流體模塊基礎模型中滲透率縱橫比通常取0.1，相對滲透率多采用VGM模型，而毛細壓力模型采用VG函數。對于力學模型，一般的數值模擬研究都假定CO2注入引起的垂向總應力不變，而地層側向變形受限，只在垂直方向上變形。

圖1 CO2咸水層封存多場耦合原理

1.2 力學本構模型

用于CO2地質封存力學模擬的巖體本構模型主要有彈性、塑性和非連續介質模型3類[20]。

(1) 彈性本構模型。不能描述復雜應力下的塑性變形狀態，但其簡單快捷，計算結果收斂，易于進行高耦合程度的多場耦合計算，因此，在CO2地質封存力學模擬中仍占較大比例[21]。

(2) 塑性本構模型。主要包括彈塑性和黏塑性本構模型。復雜應力條件下，儲層和蓋層巖體的力學特征并不遵循彈性介質的應力-應變曲線，此時選擇彈塑性本構模型能更好地刻畫其變形屈服過程，尤其當模擬儲層和蓋層的長期蠕變特征時，可使用黏塑性本構模型[22]。

(3) 非連續介質模型。針對CO2的注入可擾動初始應力場，導致新裂隙產生的問題，主要利用非連續介質模型進行儲層及蓋層的巖體力學模擬[23]。非連續介質模型能更好地刻畫CO2封存系統中儲層及蓋層的裂隙特征，因此，使用其進行裂隙啟閉以及擴展過程的模擬是未來數值模擬領域的發展趨勢。

當前的數值模擬方法主要按照以下規則分析蓋層的力學完整性：對于彈性本構模型，蓋層單元應力狀態滿足剪切或張拉破壞準則[24]；對于彈塑性本構模型，蓋層底部單元出現塑性屈服[25]；對于黏塑性本構模型，蓋層單元的塑性應變達到設定的某個臨界值[22]；對于非連續介質模型，蓋層單元節點斷開或已有裂隙發生張開或滑移[26]。

1.3 耦合分析方法

針對CO2地質封存中的多場耦合問題，存在多種數值分析方法，根據求解過程耦合程度的不同，可大致歸結為全耦合、弱耦合與單向耦合3類。

(1) 全耦合。需要同時求解熱力學、滲流力學和巖石力學這3類特征差異較大的方程，所有參數在每個迭代步計算完成后同時更新代入剛度矩陣，用于下一時間步的迭代計算。該方法結果準確，但需花費大量時間，因此，適用于尺度較小、相對簡單的幾何模型，例如僅用于熱-水動力場之間的耦合，且一般采用彈性和塑性本構模型[27]。

(2) 弱耦合。將偏微分方程拆解為流體模塊和力學模塊，每個模塊在一個時間步內先按照一定的順序序列獨自完成計算，再將結果代入對方下一時間步的計算中修正輸入參數，如此反復迭代執行直至結束。該方法計算結果精度有所降低，但算法實現相對容易，速度快，收斂性較好，已成為當前CO2地質封存研究中應用最廣泛的耦合分析方法。

(3) 單向耦合。各個模塊獨立運行，計算效率最高，收斂性最好，但力學與流體2個獨立模塊的計算結果不反饋給對方。因此，該方法多用于只關注力學而不注重滲流過程的模擬，且往往用于分析非連續介質模型[28]。

1.4 數值模擬分析軟件

能夠用于CO2地質封存流固耦合過程分析的數值模擬軟件有很多(表1)。滲流場分析可采用TOUGH2、ECLIPSE和CMG-GEM等[29]，應力場分析可采用FLAC3D、Abaqus和Petrel等。針對CO2地質封存流固耦合數值模擬研究，使用最多且效果較好的是FLAC3D和TOUGH2。FLAC3D可對CO2注入后蓋層巖石應力-應變的整個動態過程進行追蹤與反演，實現不同時空尺度的模擬及定量化研究，具有極強的可塑性與應用靈活性，同時可節約成本，兼顧各個因素的共同影響。而TOUGH2則擁有通用的多相流動方程及熱對流傳導數學模型，能夠對多相流、多組分、非等溫流體流動中的熱和水動力耦合過程進行分析，且模擬結果準確可靠。由于TOUGH2不能刻畫由水動力變化而導致的巖體變形過程，而力學過程的模擬在CO2地質封存問題研究中至關重要。因此，一些學者設計了TOUGH2-FLAC3D接口程序，對CO2地質封存中的相關問題進行了模擬研究。此外，為了對比分析不同數值模擬器的模擬結果，學者們分別采用了TOUGH2-FLAC3D、TOUGH2/ECO2N、TOUGHREACT等多種模擬器進行了研究[30]。目前，基于TOUGH2-FLAC3D的耦合分析已被成功用于模擬CO2注入咸水層后壓力的抬升、蓋層的力學性質變化等問題。

表1 數值模擬軟件對比分析

2 CO2咸水層封存中蓋層力學穩定性數值模擬研究現狀

2.1 用于巖體力學模擬分析的數值計算方法

CO2地質封存中儲層的力學響應依賴于溫度、流體壓力及應力變化，而蓋層及遠離注入點的區域力學響應則主要依賴于應力變化。注入初期，流體壓力相對較小，巖體力學響應是彈性可逆的。隨著注入壓力增大，巖體可能出現不可逆的力學損傷。當壓力足夠大時，還有可能導致剪切或拉張破壞。為模擬這一物理過程，尤其是力學方面的變化情況，需選擇合適的數值計算方法。針對地質力學問題，常用的數值計算方法有3類：適用于連續介質的(如有限單元法FEM、有限差分法FDM)、適用于非連續介質的(如離散元法DEM、非連續變形分析DDA)，以及用于連續-非連續介質的數值方法(如混合FEM/BEM)。對于這些數值方法，不能武斷地說某種方法優越于另一種。在地質力學建模中，數值方法的選擇需考慮較多因素，如研究目的、模擬的裂縫幾何體尺度等。研究顯示，基于連續介質的數值方法在CO2封存的力學建模分析中占主導地位[31]。相關研究中常將數值計算方法組合起來或鏈接到相應的模擬軟件，并使用彈性或彈塑性本構模型進行耦合分析，刻畫巖體地質力學行為。以連續介質的數值計算方法為例，THM耦合分析時序方案如圖2所示。

圖2 THM耦合分析時序方案示意圖

2.2 儲、蓋層應力變化數值模擬

針對CO2咸水層封存中儲、蓋層的力學穩定性問題，一些學者從大尺度角度，使用FLAC3D、TOUGH等模擬了實際封存工程中儲層流體壓力及上覆蓋層的應力變化情況。

(1) 隨著大量CO2的注入，儲、蓋層孔隙流體壓力均逐漸增加，壓力響應對空間離散更敏感。CO2開始注入后，壓力在注入點周圍迅速積聚，形成以注入點為中心、向四周逐漸減弱傳播的環型高壓集聚帶。壓力影響的波及速度很快，且早期壓力上升范圍快于后期[32-33]。

(2) 儲層流體壓力隨注入時間的增加而增大，隨徑向距離的增大而減小；水平和垂向的總應力都增加，且前者增量大于后者；水平和垂向有效應力都減小。CO2注入階段，儲層流體壓力初始階段變化較快，隨后逐漸穩定，且注入井附近壓降較小，當徑向距離增大到一定程度時，壓降突然增大直到降至地層初始壓力[34]。隨著CO2持續注入，多個環形壓力集聚帶逐漸相連、干擾重疊，形成一定面積的壓力抬升區[35]。

(3) 儲層壓力抬升很快傳遞至蓋層，蓋層滲透率對壓力傳播有重要影響。隨著注入量的增加，注入點附近壓力明顯大于遠處區域壓力；整個注入階段，儲、蓋層交界處有效應力變化及巖石變形最大，最易發生巖石斷裂[36]。垂向壓力的變化很大程度上依賴于蓋層的參數，如滲透系數及壓縮性[37]。從時間角度而言，低滲透蓋層的壓力積聚速度快于高滲透蓋層，注入井附近壓力的最大值隨著蓋層滲透率的增加而減小[38]。蓋層滲透率的提高雖然能夠避免儲層內部出現過高壓力，但卻增加了蓋層泄漏風險[39]。

2.3 儲、蓋層壓力增量敏感性數值模擬

對CO2注入咸水層后儲、蓋層壓力增量敏感性進行了數值模擬分析，結果如下。

(1) 以注入井處的壓力為響應變量時，孔隙度的敏感性最高[40]。

(2) 儲層滲透率縱橫比越小，注入井附近垂向壓力的增量越大；高滲透率儲層可以使壓力波動傳播得更快，且在遠離注入井的區域觀測到的壓力峰值更大[11]。

(3) 蓋層滲透率及壓縮系數的變化對儲層壓力增量有較大影響。當蓋層滲透系數變大時，垂向壓力變化較大；當蓋層厚度較小且滲透率相對較高時，局部封蓋效果不佳，滲透率和滲透率縱橫比的下降會導致壓力上升幅度增加[18]。

(4) 孔滲非均質的巖層中孔隙壓力的增量要高于均質巖層中孔隙壓力的增量，并且隨時間增加，二者之間的差異也逐漸增大；孔隙度和滲透率的非均質性會導致孔隙流體壓力大幅增加，提高CO2泄漏的風險[41]。

(5) 溫度低不利于壓力消散，當注入的CO2溫度與目標地層溫度差超過一定值時，冷卻效應可能導致注入區域出現拉張裂縫[42]。

除以上因素外，儲蓋層的幾何尺寸、巖體的本構模型及參數選擇、邊界范圍設定等也會對應力增量產生影響[43-44]。例如：邊界范圍設定引起的壓力增量隨范圍的增大而逐漸減小；定壓力邊界范圍設置過小會低估儲層壓力增量，設定過大又會增加計算量。

2.4 蓋層水力破裂風險評價數值模擬

當前，CO2地質封存中蓋層力學完整性評價主要包含2個方面：一是注入壓力擾動引起的儲集層局部高壓引發的儲、蓋層張拉破壞風險；二是復雜的地質環境及巖石物性變化導致局部應力集中，由此引發剪切破壞風險。

(1) 蓋層張拉破壞風險評價。開展CO2咸水層封存工程之前，需要測試封存場所的地應力情況，以準確控制注入壓力的上限，保障蓋層的力學完整性。通常情況下，通過開展水力壓裂等場地地應力測試和室內巖心實驗，來準確測試儲、蓋層的最小水平主應力，以此評估現場注入壓力是否會引發儲、蓋層的張拉破壞。

(2) 蓋層剪切破壞風險評價。CO2咸水層封存中，蓋層剪切破壞風險評價是在巖石力學實驗的基礎上進行的，一般采用數值模擬方法，通過研究封存系統應力場及其動態變化特征，準確模擬初始地應力場及CO2注入產生的壓力擾動，依據剪切破壞準則計算破壞指標，從而評估蓋層的剪切破壞風險。

(3) 基于地層初始應力狀態分析的蓋層水力破裂形式數值模擬。CO2封存中，初始的地應力狀態是影響巖石破壞形式的主要因素。當采用彈性本構模型時，無論地層初始應力處于拉伸還是壓縮應力狀態，發生剪切破壞的可能性都高于張拉破壞的可能性；在壓縮應力狀態下，剪切破壞面不會穿透蓋層，不影響其封存性；在拉伸應力狀態下，剪切破壞面通常會貫穿整個蓋層區域，破壞面的角度(大約為60 °)比壓縮應力狀態下(大約為30 °)更大[45]；且拉應力狀態下，蓋層底部發生拉張破壞的可能性更大。雷宏武等[11]的研究顯示，CO2的注入引起流體壓力急劇增加，地層有效應力減小，可能發生剪切破壞的方向在注入過程中變化不大，剪切破壞最有可能出現在最大速率注入點上部的蓋層，尤其是在緊靠第2注入層上部的蓋層，其次為靠近地表的位置，這與Rutqvist等[45]的分析結果相似。

3 當前研究中存在的主要問題

從總體現狀來看，流固耦合下蓋層力學穩定性研究取得了一定成果，但有關地質封存中蓋層水力破裂機理的定量表征及力學完整性評價方面的數值模擬研究還存在以下問題。

(1) 相關研究所建立的地質模型大多為概化模型，受軟、硬件條件的限制，模型中對蓋層地質特征的描述精度有待進一步提高。由于地震和測井精度一般大于50 m，而蓋層厚度通常小于50 m，因此，很難通過地震和測井資料對蓋層進行精確描述。此外，通過室內實驗和露頭得到的地質資料具有分散性、擾動性，不能完整精確表征蓋層的地質特點。由于缺少野外實際資料，目前所建數值模型對地層描述通常進行了不同程度的簡化。例如，將實際地層概化為水平地層，非均質介質概化為均質介質，且大部分模型忽略了斷層的封閉-開放性的影響、毛細壓力函數和相對滲透率模型取值變化的影響等。這些影響因素都給建模過程中對蓋層特性的精確描述帶來了很多不確定性，直接影響數值模擬結果的準確性。

(2) 由于缺乏充足的實驗數據，目前還沒有較權威的能夠準確表征CO2地質封存中蓋層巖石應力應變的本構關系。當前的數值模擬研究多采用傳統的多相流體在多孔介質中的滲流理論與控制方程，由于相對滲透率模型中參數的取值難度較大，研究者只能根據現場監測數據及為數不多的室內實驗數據推斷現有模型中各參數的取值范圍，來近似表征巖石的應力應變及多相滲流特性。雖然已有大量學者對超臨界CO2的性質進行了一定研究，但對超臨界CO2流動所具有的特殊性質方面的研究較少。且模型中的方程或定律都只適用于特定的范圍及條件，當某場地的實際地質條件與模型的適用范圍相差較大時，已有模型將不再適用。

(3) 就總體現狀而言，目前仍以搭接不同的模擬器來研究CO2地質封存中的熱流固耦合過程為主，缺乏單一的THM耦合模擬器。當前，已有的力學問題耦合效應數值模擬研究中，有的單從流場角度分析儲層封存的可行性，研究了裂縫對于滲流場的影響，而很少同時關注其對應力場的影響；有的單從力學角度考慮蓋層是否會失穩，關注了CO2注入可能誘發的裂縫產生問題，而較少關注裂縫產生后對滲流場的影響。相關研究缺乏多種機理相互影響和作用下的系統性研究，難以精確描述CO2地質封存中蓋層發生水力破裂的影響因素、裂縫擴展及CO2泄漏的真實過程，導致模擬結果與實際封存效果常常存在較大差異。

(4) 缺乏考慮地層參數不確定性的蓋層力學完整性評價模型。CO2地質封存系統中，地質參數、儲層連續性、實際工程情況等均存在較大不確定性[46]。通常情況下，儲、蓋層的物理性質、初始地應力場以及巖石力學參數是依據測井數據通過數學模型計算獲得，然而由于地質環境的復雜性以及人為因素導致的測井數據資料的不完備性，使得地層參數存在大量的不確定性。因此，模型中的參數多為特定條件下的測定結果，并不能真實反映自然條件下的復雜特征。而目前的研究中缺少地層參數不確定性對儲、蓋層力學響應影響的研究，關鍵參數取值的不確定性會直接導致數值模擬的預測結果與實際情況出現偏差。

4 未來研究趨勢

(1) 精確的地質建模對于模擬結果的正確性至關重要，為了得到更接近于實際情況的模擬結果，精確表征蓋層地質特性成為了蓋層地質建模中最基礎和最重要的問題，在后續研究中有待進一步探索，以提高對實際模型中地質特征的描述精度，為CO2地質封存提供更為精確的參考。

(2) 亟需開展能夠精確表征CO2咸水層封存中蓋層巖石應力應變的本構關系的研究，以便更加準確地描述超臨界CO2注入之后巖石物性參數的變化特征，以及注入壓力對巖石應力-應變損傷特征的影響，為蓋層水力破裂機理定量表征及密封性失效風險評價奠定良好的基礎。

(3) 有待設計單一的能夠精確表征CO2地質封存特征的THM耦合模擬器，開展多尺度條件下蓋層力學完整性數值模擬研究，定量表征水力破裂機理，基于原地應力場特征和破裂機理，獲取極限注入壓力等參數，有效預測斷裂和裂縫啟閉性，為水力破裂風險的有效評估以及蓋層力學完整性的精細描述與定量評價提供準確的數據支撐。

(4) 針對地質參數的隨機性，后續工作中開展敏感參數的不確定性分析與離散化研究是很有必要的。相關研究能夠幫助構建有效的模型，對CO2地質封存中蓋層發生水力破裂的風險概率進行全面準確預測，為蓋層力學完整性評價及CO2安全、高效封存提供指導。

5 結 論

(1) 用于CO2地質封存數值模擬研究的力學本構模型有彈性、塑型以及非連續介質模型3類。當采用彈塑性和黏塑性模型進行分析時，地層初始呈拉伸應力狀態不利于蓋層的力學完整性。

(2) 用于模擬CO2地質封存的耦合分析方法有全耦合、弱耦合以及單向耦合3類，其中，弱耦合應用最為廣泛。相應的流固耦合過程分析數值模擬軟件有TOUGH2、ECLIPSE、FLAC3D等，基于TOUGH2-FLAC3D的滲流-應力耦合分析在蓋層力學完整性研究中應用較為廣泛。

(3) 針對CO2地質封存力學問題模擬，常用的數值計算方法有適用于連續介質的、適用于非連續介質的，以及用于連續-非連續介質的數值計算方法。其中，基于連續介質的數值方法在模擬巖石破壞的過程，特別是破壞前的損傷演化分析中非常有效。

(4) CO2注入過程中，儲、蓋層巖石的孔隙流體壓力響應對空間離散較敏感；儲層流體壓力隨注入時間增加而增大，隨徑向距離增大而減小；儲層壓力抬升很快會傳遞至蓋層，儲、蓋層交界處最易發生破裂；蓋層滲透率對壓力傳播有重要影響，而孔隙度、滲透率、非均質性、模型邊界范圍設置等均對壓力增量有一定影響。

(5) 深部咸水層封存CO2是地質封存技術中最有前景的途徑，從總體現狀來看，流固耦合作用下蓋層力學完整性數值模擬研究取得了一定成果，但有關多種耦合機理和影響因素共同作用下蓋層應力應變本構關系的模擬研究、模型中邊界條件的設置與參數的準確獲取，以及模型精度的提高等都是后續研究中有待進一步探討的問題。