地質(zhì)封存二氧化碳沿?cái)鄬有孤?shù)值模擬研究

張志雄，謝 健，戚繼紅，胡立堂，張可霓，吳禮舟

(1.成都理工大學(xué)環(huán)境與土木工程學(xué)院，四川 成都 610059；2.北京師范大學(xué)水科學(xué)研究院，北京 100875； 3.勞倫斯伯克利國(guó)家實(shí)驗(yàn)室，美國(guó)加州 伯克利 94720)

利用深部咸水層進(jìn)行CO2地質(zhì)封存(GCS)被認(rèn)為是實(shí)現(xiàn)大規(guī)模減排、緩解溫室氣體效應(yīng)最具前景的技術(shù)之一[1]。隨著近年GCS工程在全球的不斷推進(jìn)和實(shí)施，其安全性和泄漏風(fēng)險(xiǎn)愈來(lái)愈引起人們重視。由于含水層不似廢棄油氣田，未有儲(chǔ)氣“經(jīng)歷”，因此其蓋層(隔水層)完整性和地層垂向組合特征是否滿足GCS密封性要求未得到充分驗(yàn)證。斷層作為影響儲(chǔ)層地質(zhì)條件復(fù)雜性的重要不確定性因素，對(duì)GCS工程的安全性、封存效果具有至關(guān)重要的影響[2]。斷層及其破碎帶往往是結(jié)構(gòu)復(fù)雜的非均質(zhì)系統(tǒng)[3]，GCS場(chǎng)地如果存在穿透主要儲(chǔ)層和蓋層的斷層，CO2極有可能沿著斷層破碎帶遷移至淺部地層甚至泄漏至地表。這不僅降低儲(chǔ)層的封存能力，而且可能導(dǎo)致上覆含水層污染、激活斷層等環(huán)境地質(zhì)效應(yīng)[4]。

研究地質(zhì)封存CO2及其儲(chǔ)層咸水在注入和封存階段沿?cái)鄬舆w移泄漏的主要特征和受控因素是極為必要的。然而，目前國(guó)內(nèi)外關(guān)于這一課題的研究還十分有限。Pruess[5～7]研究了CO2泄漏過(guò)程中的物理化學(xué)作用。鄭菲等[8]對(duì)蘇北盆地CO2地質(zhì)封存泄漏風(fēng)險(xiǎn)的全局敏感性分析，發(fā)現(xiàn)與毛細(xì)壓力有關(guān)的斷層參數(shù)對(duì)斷層中CO2總量的影響最大。董華松等[9]重點(diǎn)關(guān)注了GCS泄漏監(jiān)測(cè)技術(shù)的發(fā)展。崔振東等[4,10～11]研究了CO2泄漏的風(fēng)險(xiǎn)及后果，指出泄漏可導(dǎo)致淺部淡水資源被污染、生態(tài)系統(tǒng)和人類健康受威脅，還有可能激活斷層。總體來(lái)看，國(guó)內(nèi)學(xué)者多數(shù)停留在介紹、綜述國(guó)外研究成果和理論，缺少這一課題的專門定量研究[12]。

本研究運(yùn)用多相流達(dá)西定律，分析儲(chǔ)層流體滲流的影響因素。在此基礎(chǔ)上利用儲(chǔ)層多相流模擬軟件TOUGH2/ECO2N建立一系列場(chǎng)地泄漏模型，旨在揭示和闡明GCS流體沿?cái)鄬有孤┑闹饕獧C(jī)理或規(guī)律，評(píng)估斷層對(duì)GCS泄漏的可能影響，為場(chǎng)地評(píng)估和封存效果研究提供理論依據(jù)。

1 理論分析

CO2注入到含水層中會(huì)引起H2O-NaCl-CO2混合系統(tǒng)多種復(fù)雜的物理化學(xué)耦合作用，其中最主要的是多組分多相流體的滲流過(guò)程，這可用多相流達(dá)西定律進(jìn)行描述[13]。向含水層中注入CO2會(huì)引起孔隙壓力的抬升，這必然會(huì)改變儲(chǔ)層流體所受的壓力梯度。而系統(tǒng)毛細(xì)壓力和流體黏度的變化也會(huì)導(dǎo)致流體流動(dòng)性的變化。如果CO2所在含水層被斷層切割且斷層導(dǎo)水，則儲(chǔ)層流體在壓力梯度和密度差引起的浮力作用下必然沿?cái)鄬影l(fā)生遷移，進(jìn)而突破蓋層進(jìn)入上覆含水層。CO2與咸水在含水層乃至斷層中的遷移擴(kuò)散的過(guò)程[1,14]：

(1)

(2)

(3)

式中：t——注入時(shí)間；

Vn——流動(dòng)系統(tǒng)的體積；

M——單位體積的質(zhì)量或能量；

K——組分；

F——質(zhì)量通量；

n——表面單元dΓn的法向矢量；

q——質(zhì)量或能量源匯項(xiàng)；

φ——孔隙度；

S——飽和度；

β——相態(tài)(氣相或液相)；

ρ——密度；

X——質(zhì)量分?jǐn)?shù)；

k——巖體固有滲透率；

kr——相對(duì)滲透率；

μ——黏度；

g——重力加速度。

由式(3)可知，在相同的儲(chǔ)層環(huán)境下，增大CO2注入速率，系統(tǒng)內(nèi)氣相相對(duì)滲透率及壓力梯度Pβ均增大，這會(huì)導(dǎo)致Fβ也增大，即泄漏速率及泄漏量增大；斷層發(fā)育位置距離注入井越近，相同注入壓力條件下斷層中的壓力梯度Pβ越大，因此泄漏速率及泄漏量也增大；斷層傾角越小，泄漏途徑增長(zhǎng)，因此系統(tǒng)內(nèi)壓力梯度Pβ減小，泄漏速率及泄漏量也因此減??；斷層斷距的變動(dòng)可改變斷層兩盤的連通性，若因斷距使兩盤的連通性增大，也就相當(dāng)于k增大，因此泄漏速率及泄漏量也增大；對(duì)于斷層厚度的改變，斷層單位寬度上的Fβ不變，但厚度增大后，泄漏量會(huì)增多；增大斷層滲透性，顯然泄漏速率及泄漏量都增大。

2 模型說(shuō)明及模擬方案

2.1 地質(zhì)模型及網(wǎng)格剖分

選取一套砂巖泥巖互層的地層為模型研究場(chǎng)地，包括3層泥巖和3層砂巖(圖1)，以砂巖3為目標(biāo)儲(chǔ)層，泥巖2為主蓋層。模型水平方向以注入井為中心，徑向延伸1 000 m，地層模擬段埋深1 000～1 220 m?；A(chǔ)模型斷層發(fā)育位置距注入井300 m，斷層垂直。另外設(shè)置了一組研究?jī)A斜斷層的模型，模型水平徑向范圍為900 m，地層埋深1 000～1 200 m。所有模型均采用矩形網(wǎng)格剖分，網(wǎng)格大小0.5～50 m不等，于注入井和斷層附近加密，加密區(qū)網(wǎng)格水平間距為0.3～1.0 m。為了更好地刻畫(huà)斷層，將斷層帶分成滲透性較高的破碎強(qiáng)烈?guī)Ш蜐B透性較小的斷層破碎影響帶，而且在研究斷層滲透性的模型中還將斷層分成上下1區(qū)和2區(qū)(圖2)。

圖1 地質(zhì)概念模型及網(wǎng)格剖分圖Fig.1 Geological model and grid design

圖2 斷層帶分區(qū)刻畫(huà)示意圖Fig.2 Schematic diagram of the zonation of faults

2.2 模型方案及參數(shù)設(shè)置

為揭示可能影響GCS泄漏的斷層特征和注入方式，設(shè)計(jì)了6個(gè)方案總共15個(gè)模擬進(jìn)行研究。方案I、II、III、IV、V、VI是分別用于探討注入速率、斷層位置、斷層傾角、斷層斷距、斷層帶厚度及斷層滲透性對(duì)泄漏的影響，各模擬方案及參數(shù)設(shè)置，見(jiàn)表1。為了方便對(duì)照，表中基礎(chǔ)模型(編號(hào)1)有重復(fù)。

模型中儲(chǔ)層及蓋層的參數(shù)是參考神華CCS項(xiàng)目的主要儲(chǔ)層劉家溝組及上部泥巖層設(shè)置的，其它參數(shù)則是在儲(chǔ)層及蓋層的基礎(chǔ)上設(shè)計(jì)的；斷層的傾角設(shè)置成90°(研究?jī)A角的除外)，這主要是為了方便刻畫(huà)模型；切穿深部地層的斷裂傾角一般較大，因此在方案III中設(shè)置了大于45°的三個(gè)傾角；對(duì)于注入速率3.17 kg/s是參考神華項(xiàng)目的10×104t/a設(shè)置。概念模型中各地質(zhì)單元的水文地質(zhì)參數(shù)具體設(shè)置見(jiàn)表2。

2.3 模型邊界條件與初始條件

模型右邊側(cè)向邊界設(shè)置為常壓邊界，以保證含水層水平方向的開(kāi)放(無(wú)限含水層)，設(shè)置成常壓邊界是將邊界上的網(wǎng)格體積設(shè)置成超大體積(1051m3)實(shí)現(xiàn)的。每個(gè)模型的左側(cè)邊界分成三段，上下兩段都設(shè)置成無(wú)滲流邊界，滲透率設(shè)成很低(10-20m2)；中段設(shè)置成常壓邊界，保證CO2能正常注進(jìn)儲(chǔ)層，只需把滲透率設(shè)置成比儲(chǔ)層的滲透率略大。

模型中各地層溫度假設(shè)為42.5 ℃，模擬過(guò)程不考慮非等溫效應(yīng)。模型頂部壓強(qiáng)參考神華CCS場(chǎng)地取12.0 MPa，利用重力平衡計(jì)算得出模型的初始?jí)毫Ψ植肌｛}度也參考神華的劉家溝組咸水層取0.03%。各地層中的初始CO2質(zhì)量分?jǐn)?shù)為0。

表1 模型方案及斷層參數(shù)設(shè)置Table 1 Model schemes and associated fault parameters

注：表中V為CO2注入速率；ζ為斷層距注入井的距離；α為斷層傾角；d為斷層斷距；λ為斷層破碎帶厚度；K1x,K1y為斷層破碎影響帶水平方向滲透率；K1z為斷層破碎影響帶垂直方向滲透率；K2x,K2y為斷層破碎強(qiáng)烈?guī)椒较驖B透率；K2z為斷層破碎強(qiáng)烈?guī)Т怪狈较驖B透率。

表2 概念模型的各地質(zhì)單元物性參數(shù)設(shè)置一覽表Table 2 Physical properties of the geological units in the models

注：表中d1為主模型各地質(zhì)單元厚度；d2為輔助模型各地質(zhì)單元厚度；n為孔隙度；ρ為密度；Kx,Ky為水平方向滲透率；Kz為垂直方向滲透率；P0為進(jìn)氣毛細(xì)壓力；m為孔隙分布指數(shù)；Slr為殘余液相飽和度；Sgr為殘余氣相飽和度。

3 模型結(jié)果及分析

3.1 泄漏時(shí)間

文中以有CO2通過(guò)蓋層底部斷層位置的時(shí)刻為開(kāi)始泄漏時(shí)間，是表征CO2在儲(chǔ)層中遷移速率的一個(gè)變量。文中儲(chǔ)層及流體的性質(zhì)都一致，所以CO2從注入到泄漏的時(shí)間長(zhǎng)短，主要取決于壓力梯度Pβ的大小和斷層位置的遠(yuǎn)近。表3為6組模擬方案的泄漏時(shí)間列表，可見(jiàn)注入速率、斷層位置及斷層傾角對(duì)泄漏時(shí)間的影響最大，其次是斷層帶的滲透性，而斷層的斷距、厚度對(duì)其影響不大。注入速率由1.585 kg/s增大到6.34 kg/s，泄漏時(shí)間提前了3 853 d；斷層位置由500 m減至100 m，泄漏時(shí)間由9 821 d提前至277 d；斷層傾角由45°增至90°，泄漏時(shí)間由132 d提前至91 d，突破模型頂部時(shí)間提前了360 d；斷層滲透性放大一倍后，雖然開(kāi)始泄漏時(shí)差別不大，但突破模型頂部時(shí)間提前了3 782 d；但考慮到斷層上下部分滲透性不同時(shí)，發(fā)現(xiàn)其泄漏時(shí)間是由兩部分共同決定的，但并不是對(duì)等關(guān)系。

3.2 泄漏速率

表3 不同模擬方案CO2泄漏時(shí)間對(duì)比Table 3 Comparison of the time at which CO2 leaks for various modeling schemes

圖3 地下水及CO2的泄漏速率與時(shí)間的關(guān)系Fig.3 Leakage rate of CO2 and groundwater

泄漏速率不同于泄漏時(shí)間，它表征的是單位時(shí)間泄漏量的多少。圖3為各組模型的地下水及CO2泄漏速率曲線。需要注意的是討論泄漏速率皆為Z向的。注入速率對(duì)CO2及地下水的泄漏速率的影響顯著(圖3a)，注入速率由1.585 kg/s抬高至6.34 kg/s，CO2的泄漏速率加快了0.07～0.075 kg/s，地下水的泄漏速率加快了0.18～0.20 kg/s。斷層位置的遠(yuǎn)近對(duì)CO2泄漏速率影響較大，距注入井500 m位置的CO2泄漏速率相對(duì)于100 m的加快了0.09～0.16 kg/s；但其對(duì)地下水的影響更大，這主要是因?yàn)榈叵滤芏雀?，而且CO2驅(qū)替地下水是一個(gè)積累的過(guò)程，斷層位置越遠(yuǎn)地下水泄漏量越大。圖3c中可見(jiàn)，斷層傾角越大泄漏速率偏小，90°傾角相對(duì)于45°傾角，CO2泄漏速率小了0.18～0.34 kg/s。圖3d中的規(guī)律不夠顯著，應(yīng)該是斷距的大小會(huì)影響Z向和R向的泄漏量大小的比例。斷層帶厚度的影響很顯著，厚度越窄速率越快。斷層厚度由20 m減為10 m，地下水速率增加了0.065～0.112 kg/s，CO2泄漏速率增加了0.048～0.054 kg/s。圖3f四個(gè)曲線差別不大，斷層滲透性提升一倍，CO2泄漏速率提高了約0.005 kg/s。

3.3 CO2羽體分布及泄漏量

CO2發(fā)生泄漏后，泄漏量的大小可以用于評(píng)判CO2封存效果及封存安全性。有學(xué)者認(rèn)為，CO2的泄漏量不應(yīng)該超過(guò)注入總量的1‰[9]。圖4為六個(gè)方案10a或20a后的CO2羽體分布，水平方向?yàn)镽向，垂直方向?yàn)閆向(深度)，注入井位置在圖的左側(cè)。泄漏量的大小可大致通過(guò)CO2羽分布范圍大小來(lái)作比較，但為了更加精確的分析，圖中右側(cè)有注明儲(chǔ)層殘余CO2的比例。注入總量是注入速率與注入時(shí)間的乘積，基礎(chǔ)模型的CO2注入總量為1.99×106t。圖中可見(jiàn)CO2羽體在儲(chǔ)層中因密度小呈現(xiàn)“漏斗”狀，而在上部砂巖層因遷移先后呈現(xiàn)“倒漏斗”狀。

綜合10a或20a的CO2分布和儲(chǔ)層殘余CO2發(fā)現(xiàn)，注入速率、斷層位置對(duì)CO2泄漏量影響最大。注入速率(圖4a)為1.585kg/s時(shí)，儲(chǔ)層殘余CO2占注入總量的81.97%；速率增大一倍后，殘余CO2為60.40%；速率增大3倍后，殘余CO2為49.54%，泄漏量增多32.43%。斷層位置(圖4b)為100 m， CO2泄漏量高達(dá)注入量63.39%；而在模型5(斷層位置500 m)中還未發(fā)生泄漏。模型傾角由45°增大至90°，泄漏量相比增大了1.95%。圖4d可見(jiàn)，斷層右盤越靠下，泄漏量越大，當(dāng)斷距為40 m時(shí)，泄漏量高達(dá)注入量的43.46%。斷層帶厚度越大，泄漏量越大，厚度20 m相對(duì)于10 m的情況，泄漏量增大了5.26%(圖4e)。斷層滲透性越強(qiáng)，泄漏的越多，滲透率加倍后，泄漏量增加了2.11%(圖4f)；當(dāng)考慮到斷層上下部分滲透性不同時(shí)，發(fā)現(xiàn)其泄漏量是由兩部分共同決定的，但并不是對(duì)等關(guān)系。

圖4 不同模擬方案的CO2羽分布圖Fig.4 Comparison of the CO2 plumes for different model schemes

3.4 結(jié)果分析

綜合上述6個(gè)方案的模擬結(jié)果發(fā)現(xiàn)，研究注入速率、斷層位置的I、II兩個(gè)方案結(jié)果極其顯著。在相同的儲(chǔ)存環(huán)境下，增大CO2注入速率，系統(tǒng)內(nèi)氣相相對(duì)滲透率及壓力梯度Pβ均增大，這會(huì)導(dǎo)致泄漏速率及泄漏量增大；斷層發(fā)育位置距離注入井越近，相同注入壓力條件下斷層中的壓力梯度Pβ越大，因此泄漏速率及泄漏量也增大。

然而研究斷層傾角、斷距、厚度及滲透性的III、IV、V、VI 4個(gè)方案結(jié)果不是很顯著。斷層傾角越小，泄漏途徑增長(zhǎng)，因此系統(tǒng)內(nèi)壓力梯度Pβ減小，泄漏速率及泄漏量也因此減??；斷層斷距的變動(dòng)可改變斷層兩盤的接觸關(guān)系，極可能影響泄漏通道的連通性，若因斷距使連通性增大，也就相當(dāng)于k增大，因此泄漏速率及泄漏量也增大；對(duì)于斷層厚度的改變，系統(tǒng)內(nèi)壓力梯度Pβ差別不大，導(dǎo)致斷層單位寬度上的Fβ變化很小，但厚度增大后，泄漏量會(huì)增多；增大斷層滲透性，顯然泄漏速率及泄漏量都增大。

4 結(jié)論

(1)6組模擬方案中，研究注入速率、斷層位置的I、II 2個(gè)方案結(jié)果極其顯著。注入速率由1.585 kg/s增大至6.34 kg/s，CO2泄漏將提前3 853 d，泄漏量增多32.43%；斷層距注入井的位置由500 m減小至100 m，CO2泄漏將提前9 544 d，泄漏量增多63.39%。

(2)6組模擬方案中，研究斷層傾角、斷距、厚度及滲透性的III、IV、V、VI 4個(gè)方案結(jié)果較顯著。斷層傾角由45°增大至90°，CO2泄漏將提前41 d，泄漏量增多1.95%；斷層斷距由-20 m增大至40 m，泄漏時(shí)間差別不大，但泄漏量增多5.97%；斷層厚度由10 m增大到20 m，泄漏時(shí)間差別不大，但泄漏量增多5.26%；斷層滲透性增大一倍，泄漏時(shí)間差別不大，但泄漏量增多2.11%；斷層上下部分滲透性不同時(shí)， CO2泄漏的影響是由兩部分共同決定的，但并不是對(duì)等關(guān)系。

(3)模型結(jié)果分析表明，CO2注入方案、斷層性質(zhì)(發(fā)育位置、產(chǎn)狀、幾何形態(tài)、內(nèi)部結(jié)構(gòu))、系統(tǒng)內(nèi)巖層的組合形態(tài)對(duì)CO2泄漏均有影響。注入速率越大、斷層位置越近、斷層傾角越大、斷層越厚、滲透性越高，CO2泄漏的越早、泄漏量也越大。但斷距對(duì)CO2的泄漏主要受各巖層厚度及組合形態(tài)，與斷距大小關(guān)系不大。在這些因素中，對(duì)CO2泄漏影響最大的是注入速率和斷層位置。

[1] 謝健，張可霓，王永勝, 等.鄂爾多斯深部咸水層CO2地質(zhì)封存效果評(píng)價(jià)[J]. 巖土力學(xué), 2016, 37(1): 166-174.[XIE J, ZHANG K N, WANG Y S,etal. Performance assessment of CO2geological storage in deep saline aquifers in Ordos Basin, China[J]. Rock and Soil Mechanics, 2016, 37(1): 166-174.(in Chinese)]

[2] Chuanhe Lu, Yunwei Sun, Thomas A,etal. Uncertainty quantification of CO2leakage through a fault with multiphase and nonisothermal effects[J]. Greenhouse Gases, 2012, 2(6): 445-459.

[3] VIALLE S, DRUHAN JL, MAHER K. Multi-phase flow simulation of CO2leakage through a fractured caprock in response to mitigation strategies[J]. International Journal of Greenhouse Gas Control, 2016,44: 11-25.

[4] 崔振東, 劉大安, 曾榮樹(shù)，等. CO2地質(zhì)封存工程的潛在地質(zhì)環(huán)境災(zāi)害風(fēng)險(xiǎn)及防范措施[J].地質(zhì)論評(píng), 2011, 57(5): 700-706.[CUI Z D, LIU D A, ZENG R S,etal. Potential geological and environmental risks and its prevention measures for CO2geological storage projects[J]. Geological Review, 2011, 57(5): 700-706.(in Chinese)]

[5] PRUESS K. Integrated modeling of CO2storage and leakage scenarios including transitions between super-and subcritical conditions, and phase change between liquid and gaseous CO2[J]. Greenhouse Gases, 2011, 1(3): 237-247.

[6] PRUESS K. On CO2fluid flow and heat transfer behavior in the subsurface, following leakage from a geologic storage reservoir[J]. International Journal of Environmental Geology, 2008, 54(8): 1677-1686.

[7] PRUESS K. Leakage of CO2from geologic storage: Role of secondary accumulation at shallow depth[J]. Greenhouse Gas Control, 2008, 2(1): 37-46.

[8] 鄭菲, 施小清, 吳吉春, 等. 蘇北盆地鹽城組咸水層CO2地質(zhì)封存泄漏風(fēng)險(xiǎn)的全局敏感性分析[J]. 高校地質(zhì)學(xué)報(bào), 2012, 18(2): 232-238.[ZHENG F, SHI X Q, WU J C,etal. Global sensitivity analysis of leakage risk for CO2geological sequestration in the saline aquifer of Yancheng formation in Subei Basin[J]. Geological Journal of China Universities, 2012, 18(2): 232-238.(in Chinese)]

[9] 董華松, 黃文輝. CO2捕捉與地質(zhì)封存及泄漏監(jiān)測(cè)技術(shù)現(xiàn)狀與進(jìn)展[J]. 資源與產(chǎn)業(yè), 2010, 12(2): 123-128.[DONG H S, HUANG W H. Researach of CO2capture, geological storage and leakage technologies[J]. Resources & Industries, 2010, 12(2): 123-128.(in Chinese)]

[10] 李毅, 張可霓, 王笑雨. CO2地質(zhì)封存泄漏對(duì)淺層地下水影響的分析評(píng)價(jià)[J]. 工程勘察, 2014, 42(11): 44-50. [LI Y, ZHANG K N, WANG X Y. Assessment of the impact on shallow groundwater system by leakage of CO2geological storage[J]. Geotechnical Investigation & Surveying, 2014, 42(11): 44-50.(in Chinese)]

[11] 田地, 馬欣, 查良松, 等. 地質(zhì)封存CO2泄漏對(duì)近地表陸地生態(tài)系統(tǒng)的影響綜述[J]. 生態(tài)與農(nóng)村環(huán)境學(xué)報(bào), 2013, 29(2): 137-145. [TIAN D, MA X, ZHA L S,etal. Review of impact of CO2leakage from geologic storage on near-surface terrestrial ecological system.[J]. Journal of Ecology and Rural Environment, 2013, 29(2): 137-145. (in Chinese)]

[12] 謝健，魏寧，吳禮舟，等. CO2地質(zhì)封存泄漏研究進(jìn)展[J].巖土力學(xué), 2017,38(增刊1):181-188. [XIE J, WEI N, WU L Z,etal. Progress in leakage study of geologic CO2 storage[J]. Rock and Soil Mechanics, 2017,38(Sup1):181-188.(in Chinese)]

[13] ZHANG Keni, XU Yaqin, LING Lulu,etal. Numerical investigation for enhancing CO2injectivity in saline aquifers[J]. Energy Procedia, 2013,37: 3347-3354.

[14] Pruess K. ECO2N: A TOUGH2 Fluid Property Module for Mixtures of Water, NaCl, and CO2[R]. Berkeley: Lawrence Berkeley National Laboratory, 2005.