于恩毅 邸 元,2) 吳 輝 曹小朋 張慶福 張傳寶

* (北京大學(xué)工學(xué)院,北京 100871)

? (北京大學(xué)地球與空間科學(xué)學(xué)院,北京 100871)

** (中國(guó)石油化工股份有限公司勝利油田分公司勘探開發(fā)研究院,山東東營(yíng) 257015)

引言

二氧化碳(CO2)被認(rèn)為是主要的溫室氣體之一[1].由于大量排放,CO2在大氣中的含量不斷攀升,引起全球氣候變化,嚴(yán)重威脅人類的生存環(huán)境,減少CO2的排放已成為當(dāng)今世界關(guān)注的焦點(diǎn)問題[2-3].在可預(yù)見的未來,化石能源仍將在能源消費(fèi)體系中扮演重要角色,CO2捕集、利用和封存(CO2capture,utilization and storage,CCUS)是實(shí)現(xiàn)化石能源低碳化利用、減少CO2排放最直接和最關(guān)鍵的技術(shù)途徑[4],能夠?yàn)閷?shí)現(xiàn)碳中和目標(biāo)提供重要支撐[5].

CO2地質(zhì)封存是指通過工程技術(shù)手段將捕集的CO2注入深部地質(zhì)儲(chǔ)層,實(shí)現(xiàn)CO2與大氣長(zhǎng)期隔絕的過程[5].如圖1 所示,封存點(diǎn)一般選擇在深度800～1000 m 以下,滲透性好且上方存在低滲蓋層的儲(chǔ)層[6].儲(chǔ)層類型主要包括咸水層、油氣藏和煤層等.在所有封存類型中,深部咸水層分布廣泛,封存容量占比約為98%[5].油氣藏也是較為理想的CO2封存場(chǎng)所,因?yàn)榇嬖谠敿?xì)的地質(zhì)勘探資料和較完備的地面設(shè)施等基礎(chǔ)條件[5],并且將CO2注入油氣藏能夠利用CO2驅(qū)替達(dá)到提高原油采收率的目的[7-9].

圖1 CO2 地質(zhì)封存的途徑[10]Fig.1 Approaches of CO2 geological storage[10]

CO2在儲(chǔ)層巖體中的封存過程較為復(fù)雜,受到地層特性和地球化學(xué)反應(yīng)等影響[11],涉及溫度場(chǎng)-滲流場(chǎng)-力學(xué)場(chǎng)-化學(xué)場(chǎng)(T-H-M-C)的耦合作用[12],如圖2 所示.CO2注入地層后會(huì)導(dǎo)致較大的流體壓力積聚,導(dǎo)致有效應(yīng)力場(chǎng)發(fā)生變化,應(yīng)力場(chǎng)的改變會(huì)影響巖體孔隙度、滲透率和毛管壓力等,從而影響CO2的注入和流動(dòng)[13-14].注入的CO2為超臨界態(tài),溫度要明顯低于周圍地層溫度,二者的溫度差導(dǎo)致局部區(qū)域內(nèi)的溫度場(chǎng)發(fā)生變化,從而改變CO2流體的密度、黏度和溶解度等,影響其滲流特性,溫度變化引起的熱應(yīng)力也會(huì)直接改變巖體的應(yīng)力狀態(tài)[13].CO2易溶于水,進(jìn)而與周圍巖石礦物發(fā)生化學(xué)反應(yīng),溶解巖石或鈣化沉淀,并可能與蓋層有機(jī)質(zhì)發(fā)生反應(yīng),導(dǎo)致蓋層滲透率和孔隙度等性質(zhì)發(fā)生變化[11,15].在多孔介質(zhì)中的化學(xué)反應(yīng)速率又受到溫度、壓力、滲流速度和CO2擴(kuò)散速率等因素的影響[16].

圖2 CO2 地質(zhì)封存多場(chǎng)耦合原理[13]Fig.2 Multi-field coupled principle of CO2 geological storage[13]

CO2地質(zhì)封存的主要機(jī)制包括: (1)構(gòu)造封存,即低滲的密封蓋層限制CO2的遷移[17-18];(2)毛細(xì)封存,即CO2被孔隙結(jié)構(gòu)中的毛管力束縛[19];(3)吸附封存,即CO2吸附在黏土礦物表面[20];(4)溶解封存,即CO2與地層水或油混溶[21-22];(5) 礦化封存,即CO2與巖石發(fā)生化學(xué)反應(yīng)生成碳酸鹽礦化物[23-24].當(dāng)CO2注入到地層后,首先圈閉在油氣藏或深層鹽水層中,CO2的運(yùn)移由多相對(duì)流過程主導(dǎo),主要由注入壓力和流體密度差異引起的浮力驅(qū)動(dòng)[25].隨著時(shí)間的推移,保持自由相的CO2通過巖石中的滲流通道向上遷移,剩余的CO2移動(dòng)非常緩慢,直到最終被殘留在孔隙中、溶解在地層水中,或作為碳酸鹽巖礦物沉淀[10,26].

CO2地質(zhì)封存有其自身的風(fēng)險(xiǎn).若CO2垂向運(yùn)移至近地表區(qū)域,會(huì)造成淺層飲用水污染[27-30]、逃逸至大氣環(huán)境[31]等風(fēng)險(xiǎn).所以,在CO2注入作業(yè)前,必須對(duì)儲(chǔ)層-蓋層系統(tǒng)進(jìn)行力學(xué)分析和穩(wěn)定性評(píng)估[32-33].蓋層一般是未受干擾、低滲透、厚且橫向分布廣泛的地層,常見為頁巖、泥巖和碳酸鹽巖等[34],具有高毛細(xì)管力和突破壓力,以防止注入流體泄漏[35].在CO2封存過程中,蓋層必須承受短期的過量注入壓力和長(zhǎng)期的浮力驅(qū)動(dòng)壓力[36].斷層和裂縫帶是蓋層巖體的軟弱帶,在CO2注入過程中,蓋層中的裂縫可能發(fā)生起裂、擴(kuò)展,斷層可能被活化,發(fā)生滑移[37-38].CO2注入還可能會(huì)導(dǎo)致儲(chǔ)層中的巖體發(fā)生壓裂,裂縫可向上延伸至蓋層,甚至穿透蓋層[39].天然裂縫/斷層和壓裂裂縫形成了CO2逸出的潛在通道,對(duì)蓋層完整性構(gòu)成威脅[35,40-41].

儲(chǔ)層-蓋層系統(tǒng)多場(chǎng)耦合數(shù)值模擬是分析蓋層完整性和斷層穩(wěn)定性并進(jìn)行預(yù)警控制的關(guān)鍵核心技術(shù)[25,42].CO2地質(zhì)封存過程的數(shù)值模擬需要考慮大時(shí)間尺度和大空間尺度下的多物理過程、多場(chǎng)耦合作用、儲(chǔ)層非均質(zhì)性和流體相態(tài)變化等因素[43-45],需要研發(fā)計(jì)算效率高、計(jì)算能力強(qiáng)的CO2地質(zhì)封存數(shù)值模擬程序.本文對(duì)CO2地質(zhì)封存儲(chǔ)層-蓋層系統(tǒng)多場(chǎng)耦合數(shù)值模擬方面的研究進(jìn)行了全面綜述,介紹了CO2地質(zhì)封存過程中T-H-M-C 多場(chǎng)耦合問題的數(shù)學(xué)模型、耦合問題數(shù)值解法和封存風(fēng)險(xiǎn)分析等方面國(guó)內(nèi)外的研究進(jìn)展,對(duì)當(dāng)前我國(guó)CO2地質(zhì)封存數(shù)值模擬技術(shù)面臨的主要問題進(jìn)行了總結(jié).

1 數(shù)學(xué)模型

CO2地質(zhì)封存多場(chǎng)耦合問題的數(shù)值模擬涉及的研究對(duì)象是巖石固體骨架與孔隙流體(CO2、水、油等)組成的多相多組分系統(tǒng),其數(shù)學(xué)模型主要包括基本控制方程和相平衡計(jì)算方法.基本控制方程包括應(yīng)力平衡方程、質(zhì)量守恒方程、能量守恒方程和化學(xué)反應(yīng)等.

1.1 應(yīng)力平衡方程

根據(jù)可變形多孔介質(zhì)多相流孔隙熱彈性理論,應(yīng)力平衡方程為[12,46]

式中,ω 和 λ 為拉梅常數(shù),u為位移,P為流體壓力,Ks為巖石骨架體積模量,Km為基質(zhì)顆粒體積模量,β為骨架的體積熱膨脹系數(shù),T為溫度.本文公式中的變量均為國(guó)際單位制.

用平均應(yīng)力表述,應(yīng)力平衡方程也可寫為[47-48]

式中,v為泊松比,σm為平均應(yīng)力,α 為Biot 系數(shù).

1.2 質(zhì)量守恒方程

采用組分模型對(duì)CO2地質(zhì)封存進(jìn)行數(shù)值模擬,α組分的質(zhì)量守恒方程為

對(duì)于 α 組分,其質(zhì)量累積項(xiàng)和傳輸項(xiàng)的表達(dá)式分別為

式中,下標(biāo) β 表示相(包含油相、氣(CO2) 相和水相),? 為孔隙度,ρ 為密度,S為飽和度,表示 α 組分在 β 相中的質(zhì)量分?jǐn)?shù),Fadv為對(duì)流項(xiàng),Fdif為擴(kuò)散項(xiàng),V為質(zhì)量流速,j為擴(kuò)散速度.

巖石基質(zhì)中流體的流動(dòng)符合達(dá)西定律[49],其表達(dá)式為

式中,k為滲透率,kr為相對(duì)滲透率,μ 為黏度,g為重力加速度,D為深度.

考慮CO2擴(kuò)散,擴(kuò)散項(xiàng)可采用 Fick 定律描述

式中,D為擴(kuò)散系數(shù).

對(duì)于低速非達(dá)西滲流和高速非達(dá)西滲流,達(dá)西定律描述的滲流速度與水力梯度之間的線性關(guān)系不再成立[50].對(duì)于低滲透-致密儲(chǔ)層中的低速非達(dá)西滲流,只有驅(qū)動(dòng)壓力梯度 ?P-ρg?D大于擬啟動(dòng)壓力梯度L時(shí),流體才能流動(dòng),其表達(dá)式為[51]

對(duì)于裂縫中的高速非達(dá)西滲流,可采用Forchheimer定律描述,其表達(dá)式為

式中,χ為高速非達(dá)西滲流系數(shù).

1.3 能量守恒方程

能量守恒方程與質(zhì)量守恒方程具有相同的形式

式中,ρR為巖石密度,CR為巖石比熱,u為比內(nèi)能.

忽略輻射傳熱,熱量的流量由熱傳導(dǎo)項(xiàng)與熱對(duì)流項(xiàng)組成,其表達(dá)式為

式中,κ 為熱傳導(dǎo)系數(shù),h為比焓.

1.4 化學(xué)反應(yīng)

CO2注入地層后發(fā)生的化學(xué)反應(yīng)主要包括CO2與地層流體之間的離子交換和溶解的CO2與地層巖石之間的礦物反應(yīng)[21].對(duì)于T-H-M-C 多場(chǎng)耦合問題,化學(xué)反應(yīng)過程通過質(zhì)量和能量守恒方程中的化學(xué)反應(yīng)源項(xiàng)R來表述,其表達(dá)式為[12,26]

式中,ξ 為化學(xué)反應(yīng)計(jì)量系數(shù),r為化學(xué)反應(yīng)速率.礦物溶解和礦化反應(yīng)的化學(xué)反應(yīng)速率計(jì)算公式為[12]

式中,A為反應(yīng)接觸面積,kR為速率常數(shù),Q為化學(xué)親和性,Keq為平衡常數(shù).化學(xué)反應(yīng)速率r若為正,表示礦化沉淀;若為負(fù),表示礦物溶解.

根據(jù)不同的反應(yīng)機(jī)理,礦物溶解和沉淀不僅受到純H2O 的催化,而且受到H+和OH-的影響.因此,速率常數(shù)kR的計(jì)算公式包含3 項(xiàng)[47]

上述平衡常數(shù)和速率常數(shù)對(duì)化學(xué)反應(yīng)動(dòng)力學(xué)模型的準(zhǔn)確性至關(guān)重要,然而由于地球化學(xué)過程的復(fù)雜性,特別是相關(guān)機(jī)制的持續(xù)時(shí)間極長(zhǎng),準(zhǔn)確的測(cè)定較為困難[26].

1.5 巖石應(yīng)力-應(yīng)變本構(gòu)關(guān)系

采用雙孔隙彈性理論描述巖石的應(yīng)力-應(yīng)變關(guān)系[52]

式中,G為剪切模量,為裂縫(i=1)和基質(zhì)(i=2)中的流體和熱效應(yīng)耦合系數(shù),δij為Kronecker符號(hào).

基于等效連續(xù)介質(zhì)模型,裂縫性巖石的本構(gòu)關(guān)系為

式中,ε 為正應(yīng)變,γ 為剪應(yīng)變,σ′為有效應(yīng)力,τ 為剪應(yīng)力,GM和GF分別為基質(zhì)和裂縫的剪切模量,EM和EF分別為基質(zhì)和裂縫的楊氏模量,ECF為裂縫的壓縮模量,αM和 αF分別為巖體中基質(zhì)和裂縫的體積分?jǐn)?shù).

由于裂縫性巖石中基質(zhì)占據(jù)的體積遠(yuǎn)大于裂縫占據(jù)的體積,因此可認(rèn)為 αM≈1,式(17)中的裂縫性巖石柔度矩陣C可以寫為基質(zhì)與裂縫柔度矩陣的加權(quán)平均

式中,CM為基質(zhì)柔度矩陣,CF為裂縫柔度矩陣.

裂縫的力學(xué)性質(zhì)與巖石基質(zhì)存在很大差別,其力學(xué)參數(shù)和滲流能力會(huì)隨施加在裂縫壁面上的法向載荷產(chǎn)生動(dòng)態(tài)變化,需要通過裂縫柔度矩陣獲得裂縫壁面位移與巖石應(yīng)力間的關(guān)系,以對(duì)裂縫孔隙度、滲透率等參數(shù)進(jìn)行更新[53].對(duì)于天然裂縫在法向加載下的變形過程,Goodman 等[54]指出裂縫在法向受壓下的變形本構(gòu)關(guān)系近似為雙曲線型模型.Bandis等[55]對(duì)大量存在天然裂縫的巖塊進(jìn)行實(shí)驗(yàn),總結(jié)出裂縫受壓變形的經(jīng)驗(yàn)性雙曲線型本構(gòu)公式.

1.6 相平衡計(jì)算及輔助方程

地層條件下油-水-氣(CO2)系統(tǒng)的相平衡計(jì)算是油氣藏地質(zhì)封存CO2數(shù)值模擬研究的核心問題之一.油-水-氣(CO2)的相平衡計(jì)算主要有平衡常數(shù)法、基于Rachford-Rice 函數(shù)的閃蒸計(jì)算法和Gibbs自由能最小化法3 種計(jì)算方法[56].

平衡常數(shù)法采用基于實(shí)驗(yàn)的經(jīng)驗(yàn)公式來計(jì)算平衡常數(shù),計(jì)算速度快、穩(wěn)定性良好,雖然避免了復(fù)雜、耗時(shí)的相平衡計(jì)算,但由于缺乏嚴(yán)密的熱力學(xué)理論基礎(chǔ),計(jì)算結(jié)果往往偏差較大.基于Rachford-Rice 函數(shù)的閃蒸計(jì)算法主要用于油-氣(CO2)兩相體系的相平衡計(jì)算,同平衡常數(shù)法相比較,有較高的精度,但計(jì)算量較平衡常數(shù)法大.適用于油-水-氣(CO2) 三相體系的閃蒸計(jì)算法還有待研究和發(fā)展.Gibbs 自由能最小化方法理論嚴(yán)密、適應(yīng)性廣、計(jì)算穩(wěn)定和收斂性好,比較適合于油-水-氣(CO2)三相的相平衡計(jì)算,但是因?yàn)樾枰磸?fù)迭代,計(jì)算量比平衡常數(shù)法和閃蒸計(jì)算法都要大很多.

由于滲流場(chǎng)、溫度場(chǎng)、應(yīng)力場(chǎng)和化學(xué)場(chǎng)之間存在耦合作用,在求解T-H-M-C 多場(chǎng)耦合問題的過程中,需要實(shí)時(shí)更新計(jì)算的物性參數(shù),主要包括毛細(xì)管力、滲透率、相對(duì)滲透率、孔隙度、流體黏度和密度等.相對(duì)滲透率曲線可以采用Brooks-Corey 模型.考慮應(yīng)力場(chǎng)、溫度場(chǎng)和化學(xué)場(chǎng)對(duì)孔隙度的影響[12,53]

式中,?0為0 應(yīng)力下的初始孔隙度大小,σm為平均總應(yīng)力(壓為正),Kb為巖體的體積模量,βT為熱膨脹系數(shù),T為當(dāng)前時(shí)刻溫度,cs為固相s組分濃度,下標(biāo) 0 表示初始狀態(tài).

基質(zhì)滲透率和毛細(xì)管力可根據(jù)孔隙度變化進(jìn)行更新

式中,k0和Pc0分別為0 應(yīng)力下的初始滲透率和初始毛細(xì)管力,d為實(shí)驗(yàn)測(cè)定的參數(shù).

溫度變化會(huì)造成流體黏度的變化,一般來說隨著溫度的升高流體黏度降低,可采用經(jīng)驗(yàn)公式描述[53]

式中,A,B,C和D均為實(shí)驗(yàn)測(cè)定參數(shù).

2 耦合問題數(shù)值解法

對(duì)于CO2地質(zhì)封存的T-H-M-C 多場(chǎng)耦合問題,可采用有限元和有限體積混合的方法進(jìn)行離散,即對(duì)于式(1)的應(yīng)力平衡方程采用有限元法離散,對(duì)于式(3)質(zhì)量守恒方程和式(10)能量守恒方程采用有限體積法進(jìn)行離散.離散后得到的方程可采用全耦合、迭代耦合、弱耦合、顯式耦合和擬耦合等數(shù)值方法進(jìn)行求解,如圖3 所示.

圖3 常用的耦合問題數(shù)值解法Fig.3 Commonly used coupling solutions

全耦合方法即為全隱式的數(shù)值求解方法,該方法對(duì)離散后的應(yīng)力守恒方程、質(zhì)量守恒方程、能量守恒方程和化學(xué)反應(yīng)方程同時(shí)求解,所有參數(shù)變量需要在每個(gè)迭代步內(nèi)進(jìn)行更新[13].全耦合方法的收斂性最好,計(jì)算精度最高,但由于耦合問題的方程數(shù)量龐大,全耦合方法的計(jì)算效率往往較低.

為了提高多場(chǎng)耦合問題數(shù)值模擬的計(jì)算速度,通常采用半隱式的方法進(jìn)行求解,即優(yōu)先對(duì)一部分控制方程(一般是質(zhì)量守恒方程和能量守恒方程)進(jìn)行隱式求解,之后將求解得到的未知量作為已知量代入另一部分尚未求解的控制方程中,求解其余的未知量.顯式耦合方法、弱耦合方法和迭代耦合方法均屬于半隱式求解方法.

迭代耦合方法中,在每個(gè)時(shí)間步的牛頓迭代內(nèi)首先對(duì)質(zhì)量和能量守恒方程進(jìn)行隱式求解,將計(jì)算得到的平均孔隙壓力傳遞至應(yīng)力平衡方程中以求解巖石的變形情況,根據(jù)巖石變形程度對(duì)孔隙度、滲透率等關(guān)鍵物性參數(shù)進(jìn)行更新,開始下一迭代循環(huán),分別對(duì)質(zhì)量守恒方程、能量守恒方程和應(yīng)力平衡方程進(jìn)行計(jì)算,直至達(dá)到收斂條件進(jìn)入下一時(shí)間步.與全耦合方法相比,迭代耦合方法顯著提高計(jì)算效率[57].

在顯式耦合中,滲流場(chǎng)、溫度場(chǎng)和應(yīng)力場(chǎng)的計(jì)算獨(dú)立進(jìn)行,每個(gè)時(shí)間步中信息僅從滲流場(chǎng)、溫度場(chǎng)向應(yīng)力場(chǎng)單向傳遞一次,即在每一時(shí)間步中僅進(jìn)行一次應(yīng)力平衡方程計(jì)算和參數(shù)更新,而應(yīng)力場(chǎng)的信息則不會(huì)反饋給滲流場(chǎng)、溫度場(chǎng).顯示耦合方法是耦合程度較低的方法,計(jì)算精度有所降低,但計(jì)算效率高、收斂性好,適用于大尺度、復(fù)雜地質(zhì)條件和復(fù)雜巖體本構(gòu)模型的力學(xué)模擬[13].

弱耦合方法是在顯式耦合方法的基礎(chǔ)上進(jìn)一步減少應(yīng)力平衡方程的計(jì)算頻率,在每個(gè)時(shí)間步內(nèi)對(duì)滲溫場(chǎng)和應(yīng)力場(chǎng)進(jìn)行獨(dú)立計(jì)算,每隔數(shù)個(gè)時(shí)間步才進(jìn)行一次應(yīng)力平衡方程的計(jì)算和參數(shù)更新,如此反復(fù)循環(huán)迭代直至結(jié)束.該方法計(jì)算結(jié)果精度有所降低,但計(jì)算速度快、收斂性好、適用范圍廣,能夠耦合不同用途的程序[13].

擬耦合方法并不進(jìn)行應(yīng)力場(chǎng)求解,僅通過經(jīng)驗(yàn)公式或數(shù)據(jù)表建立孔隙度、滲透率等與孔隙壓力的關(guān)系.該方法無法對(duì)巖體的變形過程進(jìn)行模擬,因此精度有限.

3 CO2 地質(zhì)封存風(fēng)險(xiǎn)分析

分析CO2長(zhǎng)久封存過程的安全性是多場(chǎng)耦合數(shù)值模擬的關(guān)鍵應(yīng)用.圖4 展示的是在CO2地質(zhì)封存中儲(chǔ)層-蓋層可能發(fā)生的典型地質(zhì)風(fēng)險(xiǎn),主要包括CO2泄漏、地表變形、蓋層破壞和地震誘發(fā)等[39,58-61].蓋層破壞的常見形式包括斷層活化和蓋層壓裂(裂縫擴(kuò)展),斷層活化主要與斷層剪切破壞有關(guān),蓋層壓裂(裂縫擴(kuò)展)主要與蓋層拉伸破壞有關(guān)[58].

圖4 CO2 封存儲(chǔ)層-蓋層中的地質(zhì)風(fēng)險(xiǎn)示意圖[39]Fig.4 Geomechanical risks during CO2 geological storage[39]

3.1 CO2 泄漏

CO2通過蓋層泄漏的主要機(jī)理包括壓力驅(qū)動(dòng)的流體對(duì)流[62]和濃度驅(qū)動(dòng)的分子擴(kuò)散[63],如圖5 所示.斷層/裂縫和巖石基質(zhì)存在巨大的滲透率差異[64],斷層核心可視為阻礙流動(dòng)的屏障,而斷層兩側(cè)破碎的損傷區(qū)可視為沿?cái)鄬拥牧鲃?dòng)通道[35,65].在CO2的注入過程中,流動(dòng)路徑主要受裂縫網(wǎng)絡(luò)幾何形狀控制,CO2在遇到裂縫時(shí)優(yōu)先沿裂縫方向運(yùn)移,在裂縫周圍擴(kuò)散并集中分布在裂縫周圍,如圖6 所示.隨著CO2的積累,碳酸鹽巖礦物在裂縫交叉點(diǎn)沉淀,降低了裂縫網(wǎng)絡(luò)的整體連通性,低滲透帶內(nèi)的物理圈閉效應(yīng)會(huì)顯著增強(qiáng)[66].裂縫組合導(dǎo)致儲(chǔ)層滲透率各向異性,相對(duì)滲透率和毛管壓力效應(yīng)也會(huì)影響儲(chǔ)層各向異性,從而顯著影響CO2的分布[67-69].

圖5 通過斷層泄漏到地表的潛在途徑[70]Fig.5 Potential CO2 leakage path through a fault to ground surface[70]

圖6 裂縫對(duì)CO2 運(yùn)移的影響(CO2 飽和度分布)[67]Fig.6 Effect of fractures on CO2 migration (CO2 saturation distribution)[67]

CO2在多孔介質(zhì)中的擴(kuò)散速度相對(duì)于對(duì)流速度較慢.CO2的擴(kuò)散速度對(duì)流體-巖石反應(yīng)較為敏感,流體-巖石反應(yīng)可能提高(礦物溶解孔隙度增加)或降低(礦化反應(yīng)消耗CO2)擴(kuò)散速度[70].但是現(xiàn)有研究表明,擴(kuò)散并非深層儲(chǔ)層CO2的主要泄漏機(jī)制.根據(jù)Wang 等[71]的數(shù)值模擬結(jié)果,將CO2以27 MPa的注入壓力注入儲(chǔ)層(儲(chǔ)層壓力為19 MPa),如圖7所示,穩(wěn)定注入1000 年后CO2擴(kuò)散距離僅為0.93 m.此外,根據(jù)Busch 等[70]的估算,純水通過擴(kuò)散透過10 m 厚度的薄蓋層(有效擴(kuò)散系數(shù)為2.0×10-9m2/s),大約需要105年.

圖7 CO2 擴(kuò)散運(yùn)移數(shù)值模擬結(jié)果[71]Fig.7 Numerical simulation of CO2 diffusion[71]

3.2 地表變形

阿爾及利亞In Salah 二氧化碳地質(zhì)封存項(xiàng)目已經(jīng)觀測(cè)到由CO2注入導(dǎo)致的地表抬升現(xiàn)象[72-74],如圖8 所示.地表變形的分析方法主要包括解析法、半解析法和數(shù)值模擬方法[33].Fjaer 等[75]針對(duì)簡(jiǎn)化的一維橫向薄儲(chǔ)層提出計(jì)算儲(chǔ)層垂直膨脹的解析式,Xu 等[76]針對(duì)半無限空間的問題推導(dǎo)了壓力源地面變形的解析解,Selvadurai[77]提出地表巖層和深部巖體之間相互作用的基本模型,并給出地表巖層隆起的半解析式.對(duì)于實(shí)際CO2地質(zhì)封存問題中復(fù)雜巖體的變形,需要采用數(shù)值模擬方法進(jìn)行計(jì)算.

圖8 In Salah KB-502 井地面垂直位移的演化[72]Fig.8 Transient evolution of vertical ground displacement at the In Salah KB-502 well[72]

Rinaldi 等[72]通過數(shù)值模擬分析了KB-502 地表位移演化情況,并與監(jiān)測(cè)數(shù)據(jù)進(jìn)行了比對(duì),如圖9 所示.郝術(shù)仁等[78]模擬研究了不同流量和不同壓力注入CO2引起地表位移的變化規(guī)律.Siriwardane 等[45]模擬分析了蓋層斷層/裂縫帶對(duì)地表位移的影響,無斷層/裂縫時(shí)地表變形較小且對(duì)稱,有斷層/裂縫時(shí)地表位移較大且非對(duì)稱.裂縫帶滲透率較低時(shí),地面垂直位移較小;蓋層裂縫帶離注入源越近,地表變形幅度越大.

圖9 In Salah KB-502 井2006 年12 月23 日地表變形模擬結(jié)果[72]Fig.9 The simulation results of ground displacement at the In Salah KB-502 well on December 23,2006[72]

3.3 剪切破壞及斷層活化

流體壓力的演化直接影響蓋層的地質(zhì)力學(xué)穩(wěn)定性,蓋層的低滲透性阻止了超壓向蓋層快速傳播,但在靠近儲(chǔ)層界面的蓋層部分,由于流體壓力的積聚導(dǎo)致有效應(yīng)力降低,應(yīng)力狀態(tài)向失效狀態(tài)演化[79-80],如圖10 所示.根據(jù)Mohr-Coulomb 強(qiáng)度失效準(zhǔn)則(如下式),當(dāng)作用在斷層面上的剪應(yīng)力超過斷層的剪切強(qiáng)度時(shí),斷層活化[81-82]

圖10 Mohr-Coulomb 強(qiáng)度失效準(zhǔn)則[59]Fig.10 Mohr-Coulomb failure criterion[59]

表1 國(guó)外部分CO2 封存儲(chǔ)層的力學(xué)參數(shù)[59]Table 1 Stress state at several CO2 injection sites abroad[59]

CO2注入地層后流體壓力積聚、溫度場(chǎng)的變化,會(huì)改變有效應(yīng)力場(chǎng),化學(xué)反應(yīng)會(huì)降低巖石的強(qiáng)度[84].孔隙壓力的增加會(huì)導(dǎo)致莫爾圓半徑減小并左移(參見圖11 中的紅色莫爾圓).CO2到達(dá)注入井底部時(shí)溫度低于儲(chǔ)層溫度,注入井周圍的巖石發(fā)生冷卻,溫度變化誘發(fā)溫度應(yīng)力[59],應(yīng)力狀態(tài)向剪切破壞狀態(tài)轉(zhuǎn)移(參見圖11 中的藍(lán)色莫爾圓).溫度應(yīng)力主要取決于溫差、巖石基質(zhì)的體積模量和體積熱膨脹系數(shù)[85],由溫度變化 ΔT引起的孔隙壓力變化 ΔP可以根據(jù)下式進(jìn)行估算[86]

圖11 應(yīng)力狀態(tài)耦合效應(yīng)示意圖[59]Fig.11 Coupled effects on stress state[59]

式中,Δ ζ 為流體應(yīng)變,βf為流體體積熱膨脹系數(shù),為巖石基質(zhì)體積熱膨脹系數(shù),Kf為流體體積模量,Kb為巖石體積模量.

斷層除非位于注入井附近,一般不會(huì)直接受到溫度差的影響.但在注入后期,儲(chǔ)層巖石冷卻收縮,水平應(yīng)力減小,從而引起遠(yuǎn)場(chǎng)應(yīng)力的變化,導(dǎo)致遠(yuǎn)離注入井的斷層可能處于不穩(wěn)定狀態(tài)(參見圖11 中的綠色莫爾圓)[87].

3.4 拉伸破壞及裂縫擴(kuò)展

蓋層由于CO2注入發(fā)生拉伸破壞時(shí),裂縫起裂和擴(kuò)展方向依賴于最小主應(yīng)力方向和儲(chǔ)層非均質(zhì)性[88].當(dāng)孔隙壓力P超過蓋層的最小主應(yīng)力 σ3時(shí),發(fā)生拉伸破壞,裂縫開始萌生并擴(kuò)展[89-90]

當(dāng)儲(chǔ)層或蓋層巖石剪應(yīng)力超過剪切強(qiáng)度時(shí),儲(chǔ)層或蓋層已有的裂縫會(huì)發(fā)生剪切滑移[91],如圖10 所示.如果儲(chǔ)層內(nèi)部形成壓裂裂縫、裂縫不向蓋層擴(kuò)展,或者裂縫僅在蓋層有限的范圍內(nèi)延伸時(shí),則不會(huì)發(fā)生泄漏,且能夠提高CO2的注入能力[90-93].

4 數(shù)值模擬方法研究需解決的主要問題

在地下能源和資源領(lǐng)域,國(guó)際上已經(jīng)研發(fā)出了諸多用于多孔介質(zhì)多相流數(shù)值模擬的程序和軟件,能夠針對(duì)CO2地質(zhì)封存進(jìn)行T-H-M-C 多場(chǎng)耦合或部分多物理場(chǎng)耦合的模擬計(jì)算,表2 給出了目前國(guó)外主要的可用于CO2地質(zhì)封存的數(shù)值模擬器[25-26,32].針對(duì)CO2地質(zhì)封存,近期國(guó)際石油工程師協(xié)會(huì)發(fā)布了第11 個(gè)標(biāo)準(zhǔn)算例(the 11th SPE CSP)[94].該算例包括一個(gè)實(shí)驗(yàn)室尺度的二維模型、一個(gè)油藏尺度的二維和一個(gè)油藏尺度的三維地質(zhì)模型,可用于對(duì)不同的數(shù)值模擬方法和模擬器進(jìn)行對(duì)比研究.由于CO2地質(zhì)封存需要對(duì)目標(biāo)封存場(chǎng)地儲(chǔ)層-蓋層體系的三維地質(zhì)結(jié)構(gòu)進(jìn)行精細(xì)化建模,需要考慮大時(shí)間尺度和大空間尺度下的多場(chǎng)耦合作用,所以在計(jì)算 效率和計(jì)算能力方面,對(duì)數(shù)值模擬技術(shù)提出了更高的要求.除CO2地質(zhì)封存過程中復(fù)雜的相態(tài)計(jì)算和地球化學(xué)反應(yīng)涉及的本構(gòu)關(guān)系等問題以外,在數(shù)值模擬方法研究方面,還需要解決如下問題.

表2 國(guó)外部分可用于CO2 地質(zhì)封存的數(shù)值模擬器Table 2 Overview of numerical simulators for CO2 geological storage

4.1 基于角點(diǎn)網(wǎng)格的有限單元法

角點(diǎn)網(wǎng)格可以準(zhǔn)確描述地層界面、斷層面和尖滅等復(fù)雜地質(zhì)情況以及儲(chǔ)層的空間展布和屬性分布特征,因此通常都采用角點(diǎn)網(wǎng)格進(jìn)行三維地質(zhì)體精細(xì)化建模[95-96].角點(diǎn)網(wǎng)格系統(tǒng)中,通常存在重復(fù)節(jié)點(diǎn),其節(jié)點(diǎn)無法傳遞力、位移和溫度等連續(xù)性變量,相鄰網(wǎng)格之間允許錯(cuò)動(dòng)、斷開和夾層,因此基于角點(diǎn)網(wǎng)格建立的地質(zhì)模型無法直接用于三維地應(yīng)力的有限元模擬.基于角點(diǎn)網(wǎng)格數(shù)據(jù)進(jìn)行有限元網(wǎng)格的剖分,幾何處理上較為困難,且有限元網(wǎng)格數(shù)量十分巨大.因此,需要結(jié)合角點(diǎn)網(wǎng)格和有限元方法的優(yōu)點(diǎn),建立基于角點(diǎn)網(wǎng)格的地應(yīng)力有限元模擬方法,這樣既可以完整準(zhǔn)確地保留地質(zhì)模型復(fù)雜的幾何特征,也可以利用有限單元準(zhǔn)確地計(jì)算三維應(yīng)力場(chǎng)[97].

4.2 多場(chǎng)耦合的復(fù)合介質(zhì)模型

CO2注入涉及的地層區(qū)域往往包含多尺度的裂縫系統(tǒng),從長(zhǎng)度米級(jí)的小裂縫到公里級(jí)的大斷層[33].對(duì)于多尺度裂縫系統(tǒng),僅采用多重介質(zhì)模型,會(huì)喪失斷層或大裂縫準(zhǔn)確的幾何信息;若僅采用離散裂縫模型,計(jì)算量極大,不適用于實(shí)際封存問題的模擬;采用嵌入式離散裂縫模型,連續(xù)介質(zhì)雖然能夠采用結(jié)構(gòu)化網(wǎng)格,但是需要細(xì)分裂縫,并對(duì)所有不同尺度的裂縫進(jìn)行獨(dú)立表征,將導(dǎo)致網(wǎng)格數(shù)量急劇增加,造成極大的計(jì)算負(fù)擔(dān),且不適用于各向異性的連續(xù)介質(zhì).

為了更準(zhǔn)確表征裂縫型儲(chǔ)層和提高數(shù)值模擬計(jì)算效率,應(yīng)分別對(duì)不同尺度的裂縫進(jìn)行建模: 小型裂縫縫長(zhǎng)遠(yuǎn)小于網(wǎng)格尺度,對(duì)于滲流特征的影響局限在基質(zhì)網(wǎng)格內(nèi)部,可采用等效單重介質(zhì)模型建模;中等裂縫縫長(zhǎng)與網(wǎng)格尺度相當(dāng),集中分布常以裂縫簇的形式賦存,其方向性與連通性對(duì)集中分布區(qū)域的滲流特征造成較大影響,可采用雙重介質(zhì)模型建模;大型裂縫和斷層數(shù)量有限,但對(duì)于整體滲流場(chǎng)的分布具有十分重要的影響,為了準(zhǔn)確描述每一條大型裂縫或斷層的水力特征,應(yīng)采用整體嵌入式離散裂縫模型建模[98].由此形成復(fù)合介質(zhì)模型,將離散裂縫模型、多重介質(zhì)模型和整體嵌入式離散裂縫模型結(jié)合在一起,進(jìn)而在其基礎(chǔ)上建立相應(yīng)的T-H-M-C 多場(chǎng)耦合復(fù)合介質(zhì)數(shù)學(xué)模型.

4.3 分區(qū)域的有限元-有限體積混合數(shù)值離散方法

CO2地質(zhì)封存數(shù)值模擬需要考慮儲(chǔ)層的上覆蓋層、周邊地層和下部地層,建立基層-儲(chǔ)層-蓋層-地表全地層模型[87,99].CO2地質(zhì)封存場(chǎng)地尺度的地質(zhì)建模,往往達(dá)到百萬級(jí)以上的網(wǎng)格數(shù)量[100-101],模擬計(jì)算的規(guī)模十分龐大.圖12 是勝利油田G89-1 區(qū)塊的CO2地質(zhì)封存全地層模型,其網(wǎng)格數(shù)量約240 萬,待解的自由度數(shù)已達(dá)約1500 萬.

圖12 勝利油田G89-1 區(qū)塊全地層模型Fig.12 The full synthetic field model of G89-1 block in Shengli oilfield

CO2地質(zhì)封存的應(yīng)力計(jì)算區(qū)域通常遠(yuǎn)大于滲流和溫度計(jì)算區(qū)域.CO2地質(zhì)封存數(shù)值模擬中對(duì)于滲流和溫度計(jì)算具有較高的要求,儲(chǔ)層內(nèi)網(wǎng)格要保證一定的細(xì)密程度.而全地層模型應(yīng)力場(chǎng)的計(jì)算則并不需要采用相同細(xì)密程度的網(wǎng)格劃分,由于應(yīng)力平衡方程采用有限元法計(jì)算,其主變量個(gè)數(shù)也遠(yuǎn)大于滲流場(chǎng)和溫度場(chǎng)的計(jì)算,采用粗網(wǎng)格可以避免運(yùn)算資源的浪費(fèi)[54].因此,對(duì)于全地層模型的多場(chǎng)耦合問題,在滲流和傳熱計(jì)算區(qū)域內(nèi),采用細(xì)密網(wǎng)格進(jìn)行剖分(如圖13 所示),基于有限體積法對(duì)質(zhì)量守恒方程和能量守恒方程進(jìn)行數(shù)值離散;在應(yīng)力計(jì)算區(qū)域內(nèi),采用粗網(wǎng)格進(jìn)行剖分,基于有限元法對(duì)應(yīng)力平衡方程進(jìn)行數(shù)值離散.根據(jù)粗細(xì)網(wǎng)格間的映射關(guān)系,將離散方程聯(lián)立求解,從而形成分區(qū)域的有限元-有限體積混合數(shù)值離散方法.這一方法能夠在保證儲(chǔ)層滲流場(chǎng)和溫度場(chǎng)計(jì)算精度的同時(shí),降低應(yīng)力場(chǎng)求解帶來的計(jì)算負(fù)擔(dān),有效提高計(jì)算效率.

圖13 分區(qū)域網(wǎng)格劃分示意圖[54]Fig.13 Sub-regional grid division diagram[54]

4.4 高效的并行求解技術(shù)

CO2地質(zhì)封存問題涉及的時(shí)間尺度和空間尺度較大.在多場(chǎng)耦合作用的數(shù)值模擬方面,需要發(fā)展高效的并行求解技術(shù).采用并行計(jì)算技術(shù),可將一個(gè)復(fù)雜的大型計(jì)算問題分成若干相互獨(dú)立可以同時(shí)計(jì)算的子任務(wù)[102-103].常見的分解方法大體分為時(shí)間和空間上的并行,時(shí)間并行方法常采用流水線的方式實(shí)現(xiàn),空間上并行是目前絕大多數(shù)計(jì)算問題的并行解決方案,目前常采用區(qū)域分解算法[102].

對(duì)圖12 的全地層模型注CO2問題,采用并行技術(shù)進(jìn)行多場(chǎng)耦合的數(shù)值模擬,計(jì)算得到注入3.2 年時(shí)的豎向地應(yīng)力如圖14 所示.這一計(jì)算雖然采用了并行技術(shù)(Intel Xeon Platinum CPU 48 核),但模擬耗時(shí)仍達(dá)40 h.

圖14 勝利油田G89-1 區(qū)塊儲(chǔ)層豎向有效應(yīng)力模擬結(jié)果Fig.14 The simulation results of vertical effective stress at the Shengli oilfield G89-1 block

為了大幅度提高CO2地質(zhì)封存大規(guī)模模型數(shù)據(jù)處理和數(shù)值計(jì)算的效率,可基于MPI,OpenMP,CUDA 等工具研發(fā)高效的并行求解技術(shù).其中,MPI主要適用于分布式存儲(chǔ)的可拓展的并行計(jì)算機(jī)和工作站機(jī)群;OpenMP 主要適用于共享內(nèi)存多處理器系統(tǒng)和多核處理器等體系結(jié)構(gòu),且與MPI 類似,適用于粗粒度任務(wù)并行計(jì)算;而CUDA 適用于包含GPU 設(shè)備的單個(gè)計(jì)算節(jié)點(diǎn),或者與MPI/OpenMP 協(xié)同應(yīng)用于包含多個(gè)GPU 設(shè)備的集群,主要適用于細(xì)粒度數(shù)據(jù)并行計(jì)算[104].

5 結(jié)論

儲(chǔ)層-蓋層系統(tǒng)的力學(xué)分析和穩(wěn)定性評(píng)估是實(shí)現(xiàn)CO2地質(zhì)封存工業(yè)化的關(guān)鍵問題.CO2注入地層后會(huì)導(dǎo)致較大的流體壓力積聚、有效應(yīng)力場(chǎng)變化,引發(fā)CO2泄漏、地表位移、蓋層破壞和地震誘發(fā)等地質(zhì)風(fēng)險(xiǎn).蓋層破壞的常見形式包括由剪切破壞引發(fā)的斷層活化和由拉伸破壞引發(fā)的蓋層壓裂(裂縫擴(kuò)展).蓋層裂縫的起裂、擴(kuò)展和斷層的活化、滑移,形成了CO2逸出的潛在通道,對(duì)蓋層完整性構(gòu)成了重大威脅.CO2泄漏的主要機(jī)理是濃度驅(qū)動(dòng)的分子擴(kuò)散和壓力驅(qū)動(dòng)的對(duì)流,流動(dòng)路徑主要受裂縫網(wǎng)絡(luò)幾何形狀控制,但也受到地層特性和地球化學(xué)反應(yīng)等影響,涉及溫度場(chǎng)-滲流場(chǎng)-力學(xué)場(chǎng)-化學(xué)場(chǎng)的多物理場(chǎng)耦合.

CO2地質(zhì)封存數(shù)值模擬的研究對(duì)象是巖石固體骨架與孔隙流體(CO2、水和油等)組成的多相多組分系統(tǒng),其數(shù)學(xué)模型主要包括基本控制方程和相平衡計(jì)算方法等.基本控制方程包括應(yīng)力平衡方程、質(zhì)量守恒方程、能量守恒方程和化學(xué)反應(yīng)等.對(duì)離散后的控制方程耦合求解,求解方法大體可以分為全耦合、迭代耦合、弱耦合、顯式耦合和擬耦合等算法.

為了適應(yīng)CO2地質(zhì)封存儲(chǔ)層-蓋層體系精細(xì)化建模、大時(shí)間尺度和大空間尺度的多物理場(chǎng)耦合對(duì)計(jì)算模型和計(jì)算效率的要求,還需要在裂縫性介質(zhì)數(shù)值模擬方法方面進(jìn)行重點(diǎn)研究.目前,三維地質(zhì)體精細(xì)化建模通常都采用角點(diǎn)網(wǎng)格.因此,需要將角點(diǎn)網(wǎng)格和有限元方法的優(yōu)點(diǎn)相結(jié)合,建立基于角點(diǎn)網(wǎng)格的地應(yīng)力有限元模擬方法.將離散裂縫模型、多重介質(zhì)模型和整體嵌入式離散裂縫模型相結(jié)合,建立相應(yīng)的T-H-M-C 多場(chǎng)耦合復(fù)合介質(zhì)數(shù)學(xué)模型,能夠更準(zhǔn)確地表征裂縫型儲(chǔ)層并提高數(shù)值模擬的計(jì)算效率.由于CO2地質(zhì)封存全地層模型的應(yīng)力計(jì)算區(qū)域遠(yuǎn)大于滲流和溫度計(jì)算區(qū)域,采用分區(qū)域的有限元-有限體積混合數(shù)值離散方法,能夠在保證滲流場(chǎng)和溫度場(chǎng)計(jì)算精度的同時(shí),降低應(yīng)力場(chǎng)求解的計(jì)算量.由于CO2地質(zhì)封存問題涉及的時(shí)間尺度和空間尺度較大,在多物理場(chǎng)耦合方面需要發(fā)展高效的并行求解技術(shù).