基于TOUGHREACT-MP的蘇北盆地鹽城組咸水層CO2礦物封存數(shù)值模擬

莫紹星, 龍星皎, 李 瀛, 鄭 菲, 施小清, 張可霓, 趙 良

1.表生地球化學教育部重點實驗室/南京大學地球科學與工程學院, 南京 210023 2.北京師范大學水科學研究院, 北京 100875

基于TOUGHREACT-MP的蘇北盆地鹽城組咸水層CO2礦物封存數(shù)值模擬

莫紹星1, 龍星皎1, 李 瀛1, 鄭 菲1, 施小清1, 張可霓2, 趙 良1

1.表生地球化學教育部重點實驗室/南京大學地球科學與工程學院, 南京 210023 2.北京師范大學水科學研究院, 北京 100875

基于TOUGHREACT并行版，建立了蘇北盆地鹽城組下段砂巖層在碳封存中的CO2-水-巖反應二維徑向模型，并評估儲層的礦物封存潛力，探討分析了網(wǎng)格剖分精度和儲層各向異性對礦物封存模擬的影響。模擬結果表明：在CO2封存過程中，片鈉鋁石、方解石和菱鐵礦沉淀較明顯，在CO2注入5 000 a后礦物封存的比例高達34.0%；網(wǎng)格剖分精度的不同對礦物封存反應路徑?jīng)]有影響，但粗網(wǎng)格模型計算得到的礦物封存量偏高；降低儲層的kv(垂向滲透率)在短期內(nèi)會促進CO2的水平運移，有利于溶解和礦物封存；但隨著時間延長，降低kv對對流作用的抑制開始體現(xiàn)，表現(xiàn)為1 000 a后中等kv值的模型計算出的礦物封存量高于較大kv和較小kv值的模型。

CO2的礦物封存；網(wǎng)格剖分精度；各向異性；蘇北盆地

0 引言

CO2的地質(zhì)封存是減少碳排放最有效的措施之一，地下深部咸水層因其分布廣泛、人類尚不能利用而成為最有前景的CO2地質(zhì)封存的目標層[1]。目前國外已開展了廣泛的研究，取得了顯著成果[1]。雖然我國的地質(zhì)封存潛力非?？捎^[2]，但是相關研究還比較滯后，迫切需要加強。

在咸水層封存中，CO2的封存形式有構造封存、殘余氣體封存、溶解封存和礦物封存，礦物封存是最理想的封存形式[1]。由于CO2咸水層封存涉及的物理和化學過程非常緩慢，數(shù)值模擬方法是重要的研究手段，許多模擬軟件已成功運用于相關研究中[3-5]。在數(shù)值模擬中，空間離散會引起誤差，這些誤差與網(wǎng)格剖分精度有關[6]。例如當網(wǎng)格剖分太粗時，可能會高估CO2的遷移距離[7]，難以精確刻畫CO2的遷移過程[8-9]，模型短期內(nèi)會由于數(shù)值彌散而高估了CO2的溶解量[10-12]，在長期中則由于抑制了對流作用而低估CO2的溶解量[11-13]。目前對網(wǎng)格剖分精度的影響研究較少。Audigane等[13]的研究表明，提高剖分精度對模擬結果的影響較??；但由于受計算量巨大所限，他們采用的辦法是只增加注入井附近區(qū)域的剖分精度，而降低遠離注入井區(qū)域的剖分精度。

注入到儲層中的CO2在浮力作用下將向上運移，當遇到低滲透性蓋層后將向四周運移[7, 13]。在各向異性儲層中，kv(垂向滲透率)的降低，將增加CO2向上運移的阻力，CO2將不得不沿高滲透層水平運移到儲層的中部和底部。顯然，這增加了CO2與咸水的接觸時間和接觸面積，促進了CO2溶解，進而促進了礦物的溶解和新礦物沉淀[14-15]。kv降低又會抑制對流作用，主要是因為：儲層上部的CO2徑向運移距離減小，減小了對流發(fā)生的區(qū)域[16-17]；垂向滲透速率減小[16]；更多的CO2運移到儲層的中部和底部后，下部咸水的密度也會增加，使上下密度差減小[18-20]。由此可見，改變kv對CO2的水平運移和系統(tǒng)的對流作用會產(chǎn)生相反的效果(促進或抑制)，這使各向異性對封存機制的影響變得復雜。但目前在研究各向異性的影響時，大都只分析了各項異性對其中一種作用(對流作用或CO2的水平運移)的影響，并且分析的重點是對CO2溶解的影響，對礦物封存的影響研究不夠。在Zhang等[21]的研究中，當降低kv時，CO2水平運移增強，促進了溶解和礦物封存，但該研究未說明降低kv對對流作用的影響。Jahangiri等[14]、Han等[17]和Ukaegbu等[22]也發(fā)現(xiàn)kv/kh(kh指水平滲透率)值越小，CO2溶解越多，但他們均未考慮CO2的礦物封存機制。Audigane等[13]和Lindeberg等[16]的研究均表明:對流作用對CO2溶解的促進作用在系統(tǒng)的流體力學不穩(wěn)定性足夠大時才比較明顯；而kv越小，不穩(wěn)定性增加得越慢[16, 18, 23]。Kihm等[24]的研究則表明，對流作用首先促進溶解封存進而促進礦物封存。綜上所述，有必要系統(tǒng)地研究儲層各向異性對評估CO2的封存機制的影響。

蘇北盆地是含油氣沉積盆地，其深部存在體積巨大的咸水層，具有成為CO2地質(zhì)封存場所的地質(zhì)條件。本研究利用TOUGHREACT并行版軟件建立了二維徑向模型，針對蘇北盆地鹽城組下段(N1y)砂巖層內(nèi)的長期CO2-水-巖作用，分析了陽離子的反應路徑及CO2的礦物封存潛力，評估了該巖層作為碳封存儲層的可行性，并探討網(wǎng)格剖分精度和儲層的各向異性對評估礦物封存量的影響。

1 模擬軟件簡介

TOUGHREACT是基于非等溫多相流體模擬軟件TOUGH2的框架，將水流、溶質(zhì)運移和地球化學反應耦合形成的一套可變飽和地質(zhì)介質(zhì)中非等溫多相流體地球化學反應運移模擬軟件[25]，其中的ECO2N模塊已被廣泛運用于CO2地質(zhì)封存的研究中[13, 21, 25-26]。TOUGHREACT-MP是TOUGHREACT的并行版本，采用類似于TOUGH2-MP[27]的并行計算方案，引進區(qū)域分解算法對單機版進行了粗粒并行化，從而顯著提高了計算規(guī)模和效率。

2 模型的建立

2.1 區(qū)域地質(zhì)條件和模型概化

蘇北盆地內(nèi)部發(fā)育一系列北東向的南陡北緩的斷陷和凹陷，是油氣聚集的良好場所，也為CO2地質(zhì)封存提供了良好的儲存空間。鹽城組為河流環(huán)境的沉積產(chǎn)物，其下段(N1y)由3個沉積旋回組成，以粉砂質(zhì)泥巖與中粗砂巖、細礫巖互層為特征，頂部埋深一般大于1 000 m，孔隙度為23%～32%，滲透率為(0.05～1.10)×10-12m2；上段(N2y)為泥巖段，可作為區(qū)域蓋層。從埋深、巖層特征、區(qū)域蓋層封堵性、地質(zhì)圈閉條件等因素判斷，鹽城組下段是潛在的CO2封存目標層[28-29]。

在基礎模型(模型A)中，將儲層概化為徑向100 km、垂向40 m的圓柱形二維模型，儲層頂板的標高為-1 000 m。徑向上不等距遞增剖分了111層，最外一層體積設為無限大，視為一類邊界。垂向上等距剖分40層，總共剖分4 440個網(wǎng)格。此外，假設儲層上下邊界均為隔離邊界(不參與反應)。注入井位于模擬區(qū)中心，井孔位于-1 032～-1 040 m，以4 kg/s的恒定速率注入超臨界CO2，持續(xù)時間為10 a，總模擬時間為5 000 a。

為了研究網(wǎng)格剖分精度對評估CO2礦物封存量的影響，建立一個具更粗網(wǎng)格的模型(模型B)和一個具更精細網(wǎng)格的模型(模型C)，模型B在徑向和垂向上的分層數(shù)較模型A均減少一半，模型C則分別增加一倍；為了研究儲層的各向異性對評估CO2礦物封存量的影響，在模型A的基礎上另建立2個模型，分別將kv降低1個(模型D)和2個(模型E)數(shù)量級(表1)。

2.2 水文地質(zhì)參數(shù)值和儲層物質(zhì)組分的確定

表1 模型方案及研究目的

鹽城組砂巖儲層的礦物成分見文獻[32]，主要由石英、鉀長石、白云母和鈉長石組成，其次是赤鐵礦，綠泥石和綠簾石。綠泥石溶解后會產(chǎn)生碳酸鹽礦物沉淀所需的Fe2+和Mg2+，可顯著促進CO2的礦物封存[21, 24, 33-34]；綠簾石在相關研究中較少涉及，但其金屬陽離子含量較高，是碳酸鹽礦物沉淀所需二價陽離子的潛在供源。根據(jù)原生礦物的類型，參考其他學者的研究成果[13, 21, 33, 35]，也得益于TOUGHREACT并行版的運用，在模型中考慮了幾乎所有可能在碳封存過程中形成的碳酸鹽礦物和黏土礦物，次生碳酸鹽礦物包括方解石、白云石、鐵白云石、菱鎂礦、菱鐵礦和片鈉鋁石，非碳酸鹽次生礦物為高嶺石、伊利石、鈣/鈉蒙脫石和黃鐵礦(表3)。

在天然情況下，深層地下含水層中水巖相互反應已達到近似平衡狀態(tài)。為得到反應溶質(zhì)運移模型中溶液組分接近平衡的初始狀態(tài)，在進行CO2注入模擬過程之前，利用儲層的礦物組分和鹽度為2%的咸水(鹽度數(shù)據(jù)參見文獻[31]和[32])建立一個水-巖反應模型，進行1 000 a的模擬反應，將其反應1 000 a后所得的水溶液組分作為反應溶質(zhì)運移模型的初始水溶液組分。平衡后初始咸水組分見表4。

表2 模型中水文地質(zhì)參數(shù)和熱力學參數(shù)設置[30-31]

注：Slr為殘余水飽和度；Sgr為殘余氣體飽和度；p0為毛管進氣壓力；λ為指數(shù)。

表3 儲層的原生礦物、考慮的次生礦物及其反應動力學參數(shù)值[36-37]

注：假設平衡即假設該礦物和溶液可以達到局部平衡狀態(tài)。

表4 水溶液組分初始濃度和pH值

Table 4 Initial concentrations of water compositions and pH used in the simulations

組分平衡質(zhì)量摩爾濃度/(mol/kg)組分平衡質(zhì)量摩爾濃度/(mol/kg)Ca2+1.307e-5HCO-36.902e-2Mg2+2.678e-14SO2-41.000e-10Na+3.215e-1AlO2-1.131e-6K+9.461e-4Cl-2.569e-1Fe2+3.100e-7pH9.05SiO2(aq)2.509e-4

2.3 礦物反應動力學參數(shù)

在TOUGHREACT中，礦物反應(溶解或沉淀)速率所采用的模型[38-39]為

(1)

Sg為氣相CO2飽和度;b溶解相CO2為溶解相CO2質(zhì)量摩爾濃度。圖1 5 000 a后儲層中氣相CO2飽和度(a)和溶解相CO2質(zhì)量摩爾濃度(b)的分布Fig. 1 Distribution of CO2 gas saturation (a) and dissolved CO2 concentration(b)in the aquifer after 5 000 a

其中：r為礦物反應速率(mol/(L·s))；k為速率常數(shù)(mol/(m2·s))；A為每千克水的礦物比表面積(cm2/g)；K為平衡常數(shù)(無量綱)；Q為相應的離子活度積(無量綱)；θ和η為實驗測定的常數(shù)，為正值。

對于大多數(shù)礦物來說，速率常數(shù)k通常是3種機制下的綜合[33, 35]：

(2)

其中：k25為25℃時的速率常數(shù)(mol/(m2·s))；Ea為活化能(kJ/mol)；R為氣體常數(shù)(8.31 J/(mol·K))；T為絕對溫度(K)；α為組分活度(mol/L)；n為冪指數(shù)(常數(shù))；式中上、下標nu, H和OH分別代表中性、酸性和堿性機制。

表3列出了各礦物反應動力學參數(shù)的取值。其中除方解石為平衡控制外，其余礦物均為反應動力學控制。

3 模擬結果及分析

3.1 CO2飽和度的分布和pH值的變化

本文為方便起見將超臨界CO2流體稱為氣相。在超臨界CO2注入后，由于其自身的密度小于咸水密度，因而在浮力作用下向上運移，遇到蓋層阻隔后則沿其底部水平運移。圖1顯示了氣相CO2飽和度和溶解相CO2在儲層中的分布。由圖1a知，在CO2注入5 000 a后，氣相CO2羽的最大半徑約為2 600 m。在CO2羽內(nèi)，CO2溶于咸水中，使咸水密度增加，由此引起的對流作用又促進了CO2的進一步溶解(圖1b)[13, 18-19, 23]，使CO2羽本身及其附近區(qū)域形成酸性區(qū)。圖2顯示了5 000 a后儲層中pH值的分布。由圖2可見，酸性區(qū)內(nèi)的pH<5.5，酸性區(qū)外pH>10.0。

圖2 5 000 a后儲層中的pH值分布Fig. 2 Distribution of pH in the aquifer after 5 000 a

3.2 礦物的變化和陽離子反應路徑

儲層pH值的變化破壞了原有的平衡體系，導致原生礦物溶解和次生礦物沉淀，但由于溶液pH值的差異，酸性區(qū)內(nèi)外礦物的變化差別明顯，具體分析如下。

3.2.1 酸性區(qū)內(nèi)礦物的變化和陽離子反應路徑

在酸性區(qū)內(nèi)，綠簾石、綠泥石、白云母、鈉長石和赤鐵礦發(fā)生溶解。其中白云母體積分數(shù)減少最多(圖3a)，綠簾石(圖3b)、綠泥石(圖3c)和鈉長石(圖3d)溶解也較明顯，赤鐵礦溶解很少。綠簾石溶解產(chǎn)生的Fe2+和Ca2+，分別主要被菱鐵礦和方解石沉淀消耗(式(3))。綠泥石溶解產(chǎn)生Fe2+和Mg2+，F(xiàn)e2+也主要被菱鐵礦消耗(式(4))，伊利石大量沉淀消耗了Mg2+(式(4)，圖3e)，是酸性區(qū)內(nèi)沉淀最多的礦物。鈉長石溶解產(chǎn)生的Na+主要被片鈉鋁石沉淀消耗(式(5))，少部分消耗于鈉蒙脫石的沉淀。另外兩種原生礦物鉀長石和石英在酸性區(qū)內(nèi)發(fā)生沉淀，其中鉀長石的K+來自白云母(式(6))。主要反應方程式如下：

Ca2FeAl2(SiO4)3(OH)(綠簾石)+3SiO2(石英)+2K++3CO2+0.5H2O→2KAlSi3O8(鉀長石)+2CaCO3(方解石)+2H++0.25O2+FeCO3(菱鐵礦) ；

(3)

0.25Fe2.5Mg2.5Al2Si3O10(OH)8(綠泥石)+0.6KAl2(AlSi3O10)(OH)2(白云母)+0.95SiO2(石英)+0.75H++0.625CO2→K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2(伊利石)+0.625FeCO3(菱鐵礦)+0.375Mg2++0.975H2O ；

(4)

NaAlSi3O8(鈉長石)+CO2+H2O→NaAlCO3(OH)2(片鈉鋁石)+3SiO2(石英)；

(5)

KAl3Si3O10(OH)2(白云母)+ 6SiO2(石英)+2K+→ 3KAlSi3O8(鉀長石)+2H+。

(6)

體積分數(shù)變化值為負數(shù)表示溶解。圖3 5 000 a后礦物的體積分數(shù)(a--e)的變化和孔隙度Φ的分布Fig. 3 Changes of minerals in volume fraction (a-e)and distribution of porosity Φ after 5 000 a

3.2.2 酸性區(qū)外礦物的變化和陽離子反應路徑

在酸性區(qū)外，綠簾石、石英和白云母溶解較明顯，但由于pH值較高，反應強度小于酸性區(qū)內(nèi)，綠泥石則幾乎沒有溶解,鈉長石、赤鐵礦和鉀長石發(fā)生沉淀。次生礦物中，只有方解石沉淀較明顯，而在酸性區(qū)內(nèi)沉淀最多的碳酸鹽礦物片鈉鋁石則沒有沉淀，這是由于片鈉鋁石在高CO2分壓下穩(wěn)定[40]，只有在注入井附近能滿足這一條件。根據(jù)礦物的變化推知，方解石的Ca2+和赤鐵礦的Fe3+均來自綠簾石(式(7))，鉀長石的K+來自白云母(式(6))，而由于沒有礦物溶解提供Na+，所以鈉長石的沉淀消耗了溶液中原有的Na+，反應方程式為

Ca2FeAl2(SiO4)3(OH)(綠簾石)+3SiO2(石英)+2Na++2CO2+0.5H2O→2NaAlSi3O8(鈉長石)+2CaCO3(方解石)+0.5Fe2O3(赤鐵礦)+2H+。

(7)

在模擬過程中，含鎂碳酸鹽礦物鐵白云石在模擬區(qū)域內(nèi)沉淀極少，菱鎂礦和白云石則沒有沉淀，原因是伊利石作為酸性區(qū)內(nèi)沉淀最多的礦物，消耗了大量Mg2+，導致含鎂碳酸鹽礦物沉淀極少。高嶺石和鈣蒙脫石只在酸性區(qū)內(nèi)有很少量沉淀，黃鐵礦則沒有沉淀。

3.3 孔隙度和滲透率的變化

圖3f顯示5 000 a時儲層的孔隙度在模擬區(qū)域內(nèi)均降低，在酸性區(qū)內(nèi)降低幅度最大。這是由于沉淀的礦物主要是密度較小的黏土礦物伊利石、碳酸鹽礦物(片鈉鋁石、菱鐵礦和方解石)及石英等,而溶解的礦物中，重礦物綠簾石的密度較大，結果導致儲層孔隙度降低，最大降低幅度達3.2%。滲透率的變化類似于孔隙度的變化，最大降低幅度達9.3%；kv減小將減弱對流作用[16]，不利于CO2的溶解和礦物封存。

3.4 礦物封存

碳酸鹽礦物沉淀能永久固定CO2，圖4a、b顯示5 000 a時片鈉鋁石和菱鐵礦在酸性區(qū)內(nèi)沉淀均較明顯，在酸性區(qū)外沒有沉淀；圖4c、d顯示了方解石在1 000和5 000 a時體積分數(shù)的變化情況，1 000 a時方解石在酸性區(qū)內(nèi)沉淀很少，只在酸性區(qū)外沉淀較明顯，但5 000 a時在酸性區(qū)內(nèi)也大量沉淀，并且沉淀量高于酸性區(qū)外。這是由于礦物溶解緩沖了pH值，同時綠簾石不斷溶解提供豐富的Ca2+，利于方解石沉淀[36, 41]。 圖5顯示了儲層礦物封存量的空間分布及其隨時間的變化，表現(xiàn)出逐漸增加的趨勢，根據(jù)碳酸鹽礦物的沉淀情況判斷，鹽城組下段砂巖層具有較高的礦物封存潛力。儲層中所含的少量綠簾石在酸性區(qū)內(nèi)外溶解均比較明顯，產(chǎn)生的二價陽離子(Fe2+和Ca2+)主要被菱鐵礦和方解石沉淀消耗(式(3)和(7))，說明綠簾石促進了CO2的礦物封存。

a、b、d. t=5 000 a；c. t=1 000 a。圖4 片鈉鋁石(a)、菱鐵礦(b)和方解石(c-d)體積分數(shù)變化Fig. 4 Changes of dawsonite (a), siderite (b) and calcite (c-d) in volume fraction

a.t=1 000 a；b. t=5 000 a。SMCO2.儲層每立方米介質(zhì)的礦物封存量。圖5 不同時刻SMCO2的分布Fig. 5 Distribution of SMCO2 at different time

4 敏感性分析

4.1 網(wǎng)格剖分精度對評估礦物封存量的影響

數(shù)值模擬的結果中可能包含離散誤差、舍入誤差和迭代殘差[6]，本次研究擬從網(wǎng)格剖分精度的角度來研究這些誤差對評估CO2礦物封存量的影響。為此，在模型A的基礎上，建立了如表1所列的粗網(wǎng)格模型B和精細網(wǎng)格模型C。由于模型C計算量太大，因此3個模型都只模擬到1 000 a。

圖6顯示了模型B和C在1 000 a時礦物封存量的空間分布(相同時刻采用同一圖例，下同)。結合圖5a可發(fā)現(xiàn)，剖分精度越低，模擬所得的礦物封存量越大。1 000 a時，模型A、B和C礦物封存比例分別為4.5%、6.2%和2.7%。剖分最粗的模型B的礦物封存量較大的區(qū)域在水平方向上延伸約為1 900 m，但模型A和C均超過了2 000 m。這說明網(wǎng)格剖分越精細，越能精確刻畫CO2的運移，這與Juanes等[8]的研究結果類似。值得注意的是，網(wǎng)格剖分精度的不同對離子的反應路徑并未產(chǎn)生影響。

4.2 各向異性對評估礦物封存量的影響

為研究儲層各向異性對評估CO2礦物封存量的影響，在模型A的基礎上，另建立了kv/kh=0.10(模型D)和kv/kh=0.01(模型E)2個模型，每個模型的kh保持一致。圖7顯示了模型D和E在5 000 a時CO2飽和度的分布。結合圖1a可知：kv/kh越小，越多的氣相CO2運移到儲層的中部和底部；但在儲層上部，則是kv/kh值越大CO2運移越遠，這與Zhang等[21]的模擬結果類似。在3 000 a之前，kv/kh越小，溶解封存量越高(圖8)，說明水平運移對CO2溶解的促進是主要的；但3 000 a時，模型A的溶解封存量超過了模型D，說明系統(tǒng)此時的流體動力學不穩(wěn)定性已足夠大，對流作用強度增大，減小kv對對流作用強度的抑制和減小對流發(fā)生的區(qū)域的作用開始明顯體現(xiàn)。3個模型的溶解封存量在達到最大值后均緩慢減少，這與Kihm等[24]的模擬結果類似；這是由礦物封存量逐漸增加導致的。但3個模型的溶解封存量開始下降的時間不同，具中等kv值的模型D開始最早，而3 000 a后模型A下降最快，這可用對流作用對礦物封存的影響來解釋。

圖9顯示了模型D和E不同時刻的礦物封存量的空間分布，結合圖5，可看出對流作用對礦物封存的顯著影響。模型A在1 000 a以后礦物封存量較大的區(qū)域(位于酸性區(qū)及其附近區(qū)域)無論在垂向上還是水平方向，都明顯大于模型E，這種差異在5 000 a時最明顯。圖5和圖9c中礦物封存量的指狀分布也說明對流作用影響的存在，因為對流作用發(fā)生時，高密度和低密度咸水的過渡帶呈指狀分布[18-19]。礦物封存量分布的差異，是由于CO2經(jīng)過長期的運移和溶解，使系統(tǒng)產(chǎn)生足夠大的不穩(wěn)定性，對流作用已很明顯，降低kv抑制了對流作用，進而抑制了CO2和礦物的溶解，不利于礦物封存。此外，5 000 a時，在注入井周圍數(shù)百米的區(qū)域內(nèi)，模型A和D每立方米介質(zhì)的礦物封存量要比模型E多 0.5 kg左右(模型D尤其明顯)。研究表明，低pH值不利于碳酸鹽礦物沉淀[38, 41]。在模型E中，降低kv使更多CO2水平運移開后，CO2的大量溶解使溶液pH降低，雖然pH降低促進了礦物溶解，提供更豐富的陽離子，但低pH環(huán)境本身并不利于碳酸鹽礦物沉淀。由圖8可見，1 000 a后，模型D的礦物封存量超過了模型E，模型A的礦物封存量在3 000 a后快速增加，到5 000 a時超過了模型E(但仍小于模型D)。5 000 a時，模型A、D和E的礦物封存比例分別為34.0%，35.8%和33.4%。

上述結果說明，各向異性(本研究為kh固定，kv降低不同數(shù)量級)影響評估礦物封存量。在短期內(nèi)(本文為1 000 a左右)，對流作用由于系統(tǒng)的流體力學不穩(wěn)定性還未達到足夠大而并不明顯，降低kv促進氣相CO2的水平運移,從而促進了溶解和礦物封存。1 000 a后，降低kv對對流作用的的抑制開始明顯體現(xiàn)。因此在降低kv對水平運移的促進和對對流作用的抑制雙重影響下，表現(xiàn)為中等kv/kh值(模型D)相比較大kv/kh值(模型A)和較小kv/kh值(模型E)的各向異性儲層的模擬礦物封存量更高，而且在長期封存中，模型A的模擬礦物封存量超過了模型E。對流作用對礦物封存有明顯的促進作用。

圖6 1 000 a后模型B(a)和模型C(b)的SMCO2分布Fig. 6 Distribution of SMCO2 in model B (a) and C (b) after 1 000 a

圖7 5 000 a后模型D(a)和模型E(b)氣相CO2飽和度分布Fig. 7 Distribution of CO2 gas saturation in model D (a) and E (b) after 5 000 a

圖8 CO2的封存機制及其隨時間的變化Fig. 8 Time evolution of the injected CO2 in different trapping mechanisms

a、b. t=1 000 a；c、d. t=5 000 a。圖9 模型D(a, c)和模型E(b, d)的SMCO2不同時刻的分布Fig. 9 Distribution of SMCO2 in model D (a, c) and C (b, d) at different time

5 結論

本研究基于TOUGHREACT-MP對蘇北盆地鹽城組下段砂巖層CO2封存進行了二維CO2-水-巖反應數(shù)值模擬，得出了以下結論：

1) 在超臨界CO2注入后，儲層中的鈉長石、綠簾石和綠泥石溶解，促使碳酸鹽礦物片鈉鋁石、方解石和菱鐵礦沉淀，促進了CO2的礦物封存。5 000 a后礦物封存比例占34.0%，說明蘇北盆地鹽城組下段砂巖層具有較高的礦物封存潛力。由于溶解的綠簾石密度較大，而沉淀的礦物密度較小，導致儲層孔隙度和滲透率降低。

2) 對不同網(wǎng)格剖分精度的研究結果表明，剖分精度的不同對離子的反應路徑?jīng)]有影響，但粗網(wǎng)格模型計算得到的礦物封存比例偏高，數(shù)值誤差較大。因此雖然提高剖分精度會增加計算量，但為使結果盡可能精確，模型應達到一定的剖分精度。

3) 儲層的各向異性影響對CO2礦物封存量的評估。短期內(nèi)，降低kv促進了CO2的溶解和礦物封存；但隨著時間推移，降低kv對對流作用的抑制會隨著流體力學不穩(wěn)定性的增加而逐漸體現(xiàn)，進而使計算得到的礦物封存量減少。對流作用對礦物封存有明顯的促進作用。

[1] Metz B, Davidson A, de Coninck H, et al. IPCC Special Report on Carbon Dioxide Capture and Storage[R]. Geneva: IPCC, 2005.

[2] Dahowski R, Li X, Davidson C, et al. A Preliminary Cost Curve Assessment of Carbon Dioxide Capture and Storage Potential in China[J]. Energy Procedia, 2009, 1(1): 2849-2856.

[3] Krupka K M, Cantrell K J, Mcgrail B P. Ther-modynamic Data for Geochemical Modeling of Carbonate Reactions Associated with CO2Sequestration-Literature Review[R]. Richland: Pacific Northwest National Laboratory, 2010.

[4] Schnaar G, Digiulio D C. Computational Modeling of the Geologic Sequestration of Carbon Dioxide[J]. Vadose Zone Journal, 2009, 8(2): 389-403.

[5] 莫紹星, 李瀛, 龍星皎, 等. 咸水層CO2礦物封存數(shù)值模擬研究進展[J]. 地質(zhì)科技情報, 2013, 32(6): 150-158. Mo Shaoxing, Li Ying, Long Xingjiao, et al. Development of Numerical Simulation for CO2Mineral Sequestration[J]. Geological Science and Technology Information, 2013, 32(6): 150-158.

[6] 薛禹群, 謝春紅. 地下水數(shù)值模擬[M]. 北京: 科學出版社, 2007. Xue Yuqun, Xie Chunhong. Groundwater Numerical Simulation[M]. Beijing: Science Press, 2007.

[7] Beni A N, Clauser C, Erlstr?m M. System Analysis of Underground CO2Storage by Numerical Modeling for a Case Study in Malm?[J]. American Journal of Science, 2011, 311(4): 335-368.

[8] Juanes R, Spiteri E, Orr F, et al. Impact of Relative Permeability Hysteresis on Geological CO2Storage[J/OL]. Water Resources Research, 2006, 42: W12418[2013-10-23].doi:10.1029/2005WR004806.

[9] Basbug B, Gumrah F. Parametric Study of Carbon Dioxide Sequestration in Deep Saline Aquifers[J]. Energy Sources：Part A: Recovery, Utilization, and Environmental Effects, 2009, 31(3): 255-272.

[10] Green C P, Ennis-King J. Spatial Grid Correction for Short-Term Numerical Simulation Results of Carbon Dioxide Dissolution in Saline Aquifers[J]. Computational Geosciences, 2012, 16(4): 1153-1161.

[11] Ennis-King J, Gibson-Poole C M, Lang S C, et al. Long-Term Numerical Simulation of Geological Storage of CO2in the Petrel Sub-Basin, North West Australia[C]//Proceedings of the Sixth International Conference on Greenhouse Gas Control Technologies. Oxford: Pergamon Press, 2003.

[12] Ennis-King J, Lincoln P. Engineering Aspects of Geological Sequestration of Carbon Dioxide[C]//Proceedings of SPE Asia Pacific Oil and Gas Conference and Exhibition. Melbourne：Society of Petroleum Engineers, 2002.

[13] Audigane P, Gaus I, Czernichowski-Lauriol I, et al. Two-Dimensional Reactive Transport Modeling of CO2Injection in a Saline Aquifer at the Sleipner Site, North Sea[J]. American Journal of Science, 2007, 307(7): 974-1008.

[14] Jahangiri H R, Zhang D X. Effect of Spatial He-terogeneity on Plume Distribution and Dilution During CO2Sequestration[J]. International Journal of Greenhouse Gas Control, 2011, 5(2): 281-293.

[15] Anbar S, Akin S. Development of a Linear Predictive Model for Carbon Dioxide Sequestration in Deep Saline Carbonate Aquifers[J]. Computers & Geosciences, 2011, 37(11): 1802-1815.

[16] Lindeberg E, Bergmo P. The Long-Term Fate of CO2Injected into an Aquifer[J]. Greenhouse Gas Control Technologies, 2003, 1: 489-494.

[17] Han W S, Lee S Y, Lu C, et al. Effects of Permeability on CO2Trapping Mechanisms and Buoyancy-Driven CO2Migration in Saline Formations[J/OL]. Water Resources Research, 2010, 46(7): W07510[2013-10-23]. dio: 10.1029/2009wr007850.

[18] Zhang W, Li Y L, Omambia A N. Reactive Transport Modeling of Effects of Convective Mixing on Long-Term CO2Geological Storage in Deep Saline Formations[J]. International Journal of Greenhouse Gas Control, 2011, 5(2): 241-256.

[19] Pruess K, Zhang K N. Numerical Modeling Studies of the Dissolution-Diffusion-Convection Process During CO2Storage in Saline Aquifers[R]. Berkely: Lawrence Berkeley National Laboratory, 2008.

[20] Ennis-King J, Preston I, Paterson L. Onset of Con-vection in Anisotropic Porous Media Subject to a Rapid Change in Boundary Conditions[J/OL]. Physics of Fluids, 2005, 17:084107[2013-10-23]. http://dx.doi.org/10.1063/1.2033911.

[21] Zhang W, Li Y L, Xu T F, et al. Long-Term Va-riations of CO2Trapped in Different Mechanisms in Deep Saline Formations: A Case Study of the Songliao Basin, China[J]. International Journal of Greenhouse Gas Control, 2009, 3(2): 161-180.

[22] Ukaegbu C, Gundogan O, Mackay E, et al. Si-mulation of CO2Storage in a Heterogeneous Aquifer[J]. Proceedings of the Institution of Mechanical Engineers: Part A: Journal of Power and Energy, 2009, 223(3): 249-267.

[23] Xu X F, Chen S Y, Zhang D X. Convective Stability Analysis of the Long-Term Storage of Carbon Dioxide in Deep Saline Aquifers[J]. Advances in Water Resources, 2006, 29(3): 397-407.

[24] Kihm J H, Kim J M, Wang S, et al. Hydro-geochemical Numerical Simulation of Impacts of Mineralogical Compositions and Convective Fluid Flow on Trapping Mechanisms and Efficiency of Carbon Dioxide Injected into Deep Saline Sandstone Aquifers[J/OL]. Journal of Geophysical Research, 2012, 117: B06204[2013-10-23].doi:10.1029/2011JB008906.

[25] Xu T F, Sonnenthal E, Spycher N, et al. TOU-GHREACT User’s Guide: A Simulation Program for Non-Isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media, V1. 2.1[R]. Berkely: Lawrence Berkeley National Laboratory, 2008.

[26] 許天福, 金光榮, 岳高凡, 等. 地下多組分反應溶質(zhì)運移數(shù)值模擬:地質(zhì)資源和環(huán)境研究的新方法[J]. 吉林大學學報: 地球科學版, 2012, 42(5): 1410-1425. Xu Tianfu, Jin Guangrong, Yue Gaofan, et al. Surface Reactive Transport Modeling: A New Research Approach for Geo-Resources and Environments[J]. Journal of Jilin University: Earth Science Edition, 2012, 42(5): 1410-1425.

[27] Zhang K N, Wu Y S, Pruess K. User’s Guide for TOUGH2-MP-A Massively Parallel Version of the TOUGH2 Code[R]. Berkely: Lawrence Berkeley National Laboratory, 2008.

[28] 祝厚勤, 朱煜, 鄭開富. 蘇北盆地鹽城組天然氣藏成藏條件及控制因素探討[J]. 海洋地質(zhì)動態(tài), 2003,19(9): 22-26. Zhu Houqin, Zhu Yu, Zheng Kaifu. Pooling Conditions and Controlling Factors of Gas Pools in Yancheng Formation of Subei Basin[J]. Marine Geology Frontiers, 2003, 19(9): 22-26.

[29] 翟光明, 王慎言, 邱中建, 等. 中國石油地質(zhì)志: 卷八[M]. 北京: 石油工業(yè)出版社, 1992. Zhai Guangming, Wang Shenyan, Qiu Zhongjian, et al. Petroleum Geology of China: Vol 8[M]. Beijing: Petroleum Industry Press, 1992.

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

[31] 鄭菲, 施小清, 吳吉春, 等. 蘇北盆地鹽城組咸水層CO2地質(zhì)封存泄漏風險的全局敏感性分析[J]. 高校地質(zhì)學報, 2012, 18(2): 232-238. Zheng Fei, Shi Xiaoqing, Wu Jichun. et al. 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.

[32] 江蘇省地質(zhì)局石油普查隊蘇北專題隊. 蘇北中新生界巖礦特征及其在地層對比上的意見[R]. 南京: 江蘇省地質(zhì)資料館, 1961. Subei Oil Survey Team in Geological Bureau of Jiangsu Province. The Petrological and Mineralogical Characteristics of the Mesozoic and Cenozotic Stratum and Views on the Stratigraphic Correlation in North Jiangsu Province[R]. Nanjing: Geological archive of Jiangsu Province, 1961.

[33] Gundogan O, Mackay E, Todd A. Comparison of Numerical Codes for Geochemical Modelling of CO2Storage in Target Sandstone Reservoirs[J]. Chemical Engineering Research and Design, 2011, 89(9): 1805-1816.

[34] Xu T F, Apps J A, Pruess K, et al. Numerical Modeling of Injection and Mineral Trapping of CO2with H2S and SO2in a Sandstone Formation[J]. Chemical Geology, 2007, 242(3): 319-346.

[35] Xu T F, Apps J A, Pruess K. Mineral Sequestration of Carbon Dioxide in a Sandstone-Shale System[J]. Chemical Geology, 2005, 217(3): 295-318.

[36] Lasaga A C, Soler J M, Ganor J, et al. Chemical Weathering Rate Laws and Global Geochemical Cycles[J]. Geochimica et Cosmochimica Acta, 1994, 58(10): 2361-2386.

[37] Steefel C I, Lasaga A C. A Coupled Model for Tran-sport of Multiple Chemical Species and Kinetic Precipitation/Dissolution Reactions with Application to Reactive Flow in Single Phase Hydrothermal Systems[J]. American Journal of Science, 1994, 294(5): 529-592.

[38] Palandri J L, Kharaka Y K. A Compilation of Rate Parameters of Water-Mineral Interaction Kinetics for Application to Geochemical Modeling[R]. Denver: US Geological Survey, 2004.

[39] Xu T F, Sonnenthal E, Spycher N, et al. TOU-GHREACT：A Simulation Program for Non-Isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media: Applications to Geothermal Injectivity and CO2Geological Sequestration[J]. Computers & Geosciences, 2006, 32(2): 145-165.

[40] Worden R H. Dawsonite Cement in the Triassic Lam Formation, Shabwa Basin, Yemen: A Natural Analogue for a Potential Mineral Product of Subsurface CO2Storage for Greenhouse Gas Reduction[J]. Marine and Petroleum Geology, 2006, 23(1): 61-77.

[41] Druckenmiller M L, Maroto-Valer M M. Carbon Se-questration Using Brine of Adjusted pH to Form Mineral Carbonates[J]. Fuel Processing Technology, 2005, 86(14): 1599-1614.

Numerical Modeling of CO2Sequestration in the Saline Aquifer of Yancheng Formation in Subei Basin Using TOUGHREACT-MP

Mo Shaoxing1,Long Xingjiao1,Li Ying1,Zheng Fei1,Shi Xiaoqing1,Zhang Keni2,Zhao Liang1

1.KeyLaboratoryofSurficialGeochemistry,MinistryofEducation/SchoolofEarthSciencesandEngineering,NanjingUniversity,Nanjing210023,China2.CollegeofWaterSciences,BeijingNormalUniversity,Beijing100875,China

The two-dimensional radial model of CO2-water-rock reaction in carbon sequestration of Yancheng Group under section sandstone layer in Subei basin is built based on TOUGHREACT parallel version. The mineral sequestration potential in reservoir is evaluated. The influence of grid subdivision precision and anisotropy of reservoir on mineral sequestration simulations is analyzed. The simulation results show that dawsonite, calcite and siderite, significantly precipitated in the process of CO2sequestration. The total mass of CO2sequestrated by mineral trapping is as high as 34.0% after 5 000 a. The grid resolution has little impact on the reaction path of minerals sequestration, however, the total amount of CO2mineral sequestration with coarser grid is overestimated comparing to that with finer grid. The discretion of the vertical permeability (kv) would lead to a more rapid horizontal migration of injected CO2to the middle and bottom, which results in an incretion in solubility trapping and mineral trapping for a short period of time. However, for a long term, with the increasing instability of the system due to CO2dissolution, the vertical convective mixing becomes significant, it turns out to be a significant inhibition of convection by decreasingkv. It shows that the mineral sequestration calculated by middlekvvalue model is higher than the largerkvand smallerkvvalue model after 1 000 a.

mineral sequestration of CO2; grid subdivision precision; anisotropy; Subei basin

10.13278/j.cnki.jjuese.201405208.

2013-12-23

江蘇省自然科學基金項目(BK2012313)；國家大學生創(chuàng)新訓練計劃項目(G1210284016)；江蘇省產(chǎn)學研聯(lián)合創(chuàng)新資金計劃項目(BY2010136)

莫紹星(1992--)，男，布依族，主要從事CO2地質(zhì)封存技術研究，E-mail:njujinchun@163.com

施小清(1979--)，男，副教授，主要從事反應溶質(zhì)運移模擬研究，E-mail:shixq@nju.edu.cn。

10.13278/j.cnki.jjuese.201405208

P641.69

A

莫紹星, 龍星皎, 李瀛, 等. 基于TOUGHREACT-MP的蘇北盆地鹽城組咸水層CO2礦物封存數(shù)值模擬.吉林大學學報:地球科學版,2014,44(5):1647-1658.

Mo Shaoxing, Long Xingjiao, Li Ying, et al. Numerical Modeling of CO2Sequestration in the Saline Aquifer of Yancheng Formation in Subei Basin Using TOUGHREACT-MP.Journal of Jilin University:Earth Science Edition,2014,44(5):1647-1658.doi:10.13278/j.cnki.jjuese.201405208.