郭朝斌，張可霓，魯維豐，凌璐璐，郭振東

（1.北京師范大學水科學研究院，北京 100875；2.神華鄂爾多斯煤制油分公司，內(nèi)蒙古 鄂爾多斯伊旗 017209）

溫室氣體排放增加導致溫室效應加劇，目前全球變暖已經(jīng)成為國際社會關(guān)注的熱點問題。CO2的封存問題就是在這樣的背景下提出來的，CO2的捕獲和封存技術(shù)目前已經(jīng)成為一種國際公認的減排措施，被認為有巨大的減排潛力［1］。地下咸水層在世界各地廣泛分布，而且規(guī)模大，是切實可行和最具發(fā)展前景的減排途徑，被認為是長期地質(zhì)封存的最有效方法之一［2］。

近年來國內(nèi)外很多學者對CO2地質(zhì)封存的數(shù)值模擬陸續(xù)展開了研究。R.Juanes等［3］通過對CO2捕獲機制和相對滲透率的研究及數(shù)值模擬，表明交替注入水和CO2可以提高注入效率。注入的水使大面積CO2暈分裂，增大吸濕過程的時間和面積從而增強CO2的捕獲和固定，但這種方法會導致井孔底部壓力增大，受到蓋層的封閉性能、管理和經(jīng)濟方面等的限制。理解大體積CO2暈的流體力學行為對提高咸水層中CO2地質(zhì)封存技術(shù)非常重要，Christine Doughty［4］對美國加州圣華金河峽谷南部Kimberlina發(fā)電廠進行數(shù)學建模研究，利用TOUGH2對傾斜地層中CO2地質(zhì)封存大規(guī)模初步試驗中CO2暈的演變進行模擬，結(jié)果表明在4年注入周期后，注入的100×104t CO2中約20%為溶液相，25年后，CO2暈達到穩(wěn)定，溶液相態(tài)比例上升到38%，62%為超臨界狀態(tài)，其中只有3%為自由態(tài)。隨著時間繼續(xù)，暈的演變僅表現(xiàn)為緩慢的溶解作用，并伴有飽和溶解CO2咸水的對流。國內(nèi)對提高CO2封存能力、安全評估等方面研究中，許雅琴等［5］通過模擬提高低滲地層中滲透率等探討如何提高咸水層CO2封存注入率，但對滯后效應在CO2封存方面的研究相對較少。

1 方法和理論

1.1 超臨界狀態(tài)CO2的吸濕和疏干過程

將超臨界狀態(tài)CO2封存到深部咸水層會涉及到諸多多相流問題，注入的CO2為非潤濕相，水及少量溶解到水中的CO2作為潤濕相。由于潤濕相和非潤濕相流體性質(zhì)不同，運移過程中因接觸角不同和殘余氣相的產(chǎn)生會對流體運移產(chǎn)生影響。注入時CO2暈連續(xù)增長，超臨界狀態(tài)CO2(表征為氣相)作為非潤濕相驅(qū)替含水層中潤濕相(咸水)，即疏干過程。對于隨后停止注入階段，當CO2暈在浮力作用下向上運移(圖1)，不同位置經(jīng)歷過程不同，在運移前端為疏干過程，在暈尾端，潤濕相驅(qū)替非潤濕相，為吸濕過程，有部分殘余CO2被捕獲，變?yōu)椴豢蛇\動的殘余氣相，進而影響相對滲透率和毛細壓力的計算。不同位置和過程中殘余氣相不同，使得相同飽和度下相對滲透率和毛細壓力不相同，在特征曲線上表現(xiàn)出滯后現(xiàn)象。

圖1 CO2向上運移過程暈后邊緣殘余示意圖［3］Fig.1 Schematic diagram showing the trail of residual CO2as the plume migrates upward

1.2 數(shù)值模擬方法

數(shù)學模型和數(shù)值模擬在評估地質(zhì)封存CO2可行性上具有重要的作用，它們是設(shè)計和實施CO2地質(zhì)處置必須的工具［6］?？梢栽u估預測地層的封存能力，評估CO2地質(zhì)儲存的可行性和可靠性，解釋CO2在地層中的運移變化行為以及在注入到地層中后對其運移過程進行監(jiān)測和分析。

在大多數(shù)通用多相流數(shù)值模擬軟件中，相對滲透率和毛細壓力的計算函數(shù)僅與當前飽和度有關(guān)，然而實驗表明［7］，相對滲透率和毛細壓力不僅和當前區(qū)域的飽和度有關(guān)，還和流體經(jīng)歷的過程及歷史飽和度有關(guān)。盡管滯后效應特征曲線在石油工業(yè)已經(jīng)被應用，但是在包氣帶、地熱儲藏工程、核廢料地質(zhì)處理以及CO2地質(zhì)封存等這些兩相流的問題中很少被應用到［8～9］。

目前在CO2地質(zhì)封存方面，應用較為普遍的模擬軟件為TOUGH2及2008年正式公開發(fā)布用于大規(guī)模并行計算的并行版本TOUGH2-MP［10］。目前成功應用最大的模型是東京灣的CO2地質(zhì)處置模型，其單元數(shù)超過1 千萬個［11］。TOUGH+ 是繼 TOUGH2［12］之后的新一代TOUGH家族的模擬器，在繼承了TOUGH2所有功能外，主要在數(shù)組處理效率方面提高模擬的計算效率，同時也提高了一些模擬功能。TOUGH+CO2是TOUGH+系列中模擬CO2封存的模塊［13］，最新版本TOUGH+CO2中，在 iTOUGH2［14］中關(guān)于滯后效應模塊的基礎(chǔ)上改進更新了滯后效應模擬的功能。

1.3 TOUGH+CO2中相對滲透率相關(guān)理論

滯后效應特征曲線不僅取決于當前飽和度，還和歷史飽和度以及經(jīng)歷的過程即疏干過程(非潤濕相替代潤濕相)或吸濕過程(潤濕相替換非潤濕相)相關(guān)［15～17］。

CO2深部咸水層封存中主要涉及多相流系統(tǒng)為H2O－CO2－NaCl系統(tǒng)，超臨界狀態(tài)CO2在TOUGH+CO2中表征為非水相，即氣相，溶解于水中的CO2表征為溶液相；計算系統(tǒng)中流體相的相對滲透率的常用方法為 van Genuchten［18］公式。

考慮相對滲透率滯后現(xiàn)象時，TOUGH+CO2采用Parker和Lenhard根據(jù)無滯后效應van Genuchten方程修改得到的公式［14］，

式中:krl——液相相對滲透率，無量綱；

krg——氣相相對滲透率，無量綱；

式中:Sl——液相飽和度，無量綱；

Slr——液相殘余飽和度，無量綱；

Sls——飽和水飽和度，無量綱；

Sgr——殘余氣相飽和度，無量綱。

圖2 歷史飽和度對相對滲透率的影響Fig.2 Effect of the history saturation on relative permeability

1.4 TOUGH+CO2中毛細壓力相關(guān)理論

在計算H2O－CO2系統(tǒng)中流體相的毛細壓力時經(jīng)常用到的方法為 van Genuchten［18］函數(shù)。

式中:Pcap——毛細壓力(Pa)；

P0——進入毛細壓力(Pa)。

考慮毛細壓力滯后現(xiàn)象時采用基于van Genuchten方程修改得到的公式［14］，

式中:p——流體流動過程，疏干過程(d)或吸濕過程(w)；

mp——函數(shù)適配參數(shù) mp=(nγ－1)/nγ，n、γ 為適配參數(shù)。

2 應用

2.1 模型基本參數(shù)

通過對案例模型的模擬，可以分析相對滲透率和毛細壓力滯后效應對CO2咸水層地質(zhì)封存的影響，為了研究方便，假設(shè)模型為理想條件模型。

進行CO2深部咸水層封存模擬時，需要給出的初始條件有壓力、溫度、鹽度以及流體中CO2的質(zhì)量分數(shù)。模型模擬的咸水層在地下1200m深，模擬范圍為1km×1km，模型厚度為100m，模型上下邊界為低滲透地層，其中蓋層為5m厚泥巖，注入層為85m厚砂巖，下伏層為10m厚泥巖(圖3)。模型共剖分成64000個網(wǎng)格，其中注入井附近網(wǎng)格進行加密處理。四周邊界條件視為無流量邊界，對于封閉系統(tǒng)，即儲層被完全封閉，其CO2儲存量由在儲層中壓力積聚引起孔隙體積的膨脹和咸水密度的增加來得到。

圖3 模型范圍XZ剖面示意圖Fig.3 Profile of the model in the XZ direction

此模型中只考慮等溫過程，初始溫度為37℃，模型的頂部邊界壓力為1.20×107Pa，初始壓力分布根據(jù)靜水壓力平衡計算得到。超臨界狀態(tài)下CO2的性質(zhì)根據(jù)實驗數(shù)據(jù)加入到程序，其他關(guān)鍵物理參數(shù)見表1。

表1 模型采用關(guān)鍵物理參數(shù)Table 1 Key properties of rocks in the model

方案1:不考慮滯后現(xiàn)象對CO2封存的影響，在計算中，相對滲透率和毛細壓力的計算方式均選擇van Genuchten函數(shù)。模型區(qū)域中介質(zhì)的殘余氣相飽和度Sgr為0.177，整個含水層中的殘余氣相飽和度均為固定Sgr。在數(shù)值模擬中，根據(jù)固定值Sgr，計算相對滲透率特征曲線，在驅(qū)替過程中(包括疏干和吸濕過程)，相對滲透率均沿此曲線變化，同一流體飽和度下對應相同相對滲透率。毛細壓力的計算亦是如此。

方案2:考慮滯后現(xiàn)象對CO2封存的影響。相對滲透率和毛細壓力計算函數(shù)選擇修改后的公式。模型區(qū)域中殘余氣相飽和度的最大可能值Sgrmax為0.177。在數(shù)值計算中，根據(jù)驅(qū)替過程轉(zhuǎn)換點的氣相飽和度和經(jīng)歷的過程計算相應的殘余氣相飽和度，然后通過插值得到掃描曲線，與主線不同。殘余液相飽和度Slr為0.2，krgmax為 1。

CO2注入以0.926 kg/s(約3×104t/a)恒定速率注入4年，為觀察停止注入后的運移狀態(tài)，在第5～50年停止注入，模擬時間總共為50年。

2.2 結(jié)果分析

2.2.1 對封存量的影響

在模擬時間為50年時，99%的CO2存儲于注入層中，只有1%存在蓋層和底部地層中(表2)。因浮力作用CO2向上運移，由于蓋層滲透率低，封閉性好，CO2大部分被封存在注入層中。所以兩個方案中封存的CO2總量基本一致，但CO2存在相態(tài)比例不同。

表2 模擬時間為50年時各層CO2含量Table 2 Mass of CO2in various forms at 50 years (kg)

2.2.2 注入CO2的存在形式及其分布

圖4～5為兩種方案中封存CO2相態(tài)比例變化。氣相是指整個地層中超臨界狀態(tài)CO2，包含可移動的自由態(tài)CO2和捕獲封存不可移動氣相CO2。隨著模擬時間的增加，CO2和咸水層接觸時間增加，氣相CO2逐漸溶解于溶液中成為溶液相。

圖4 無滯后效應模型整個地層中各相態(tài)CO2質(zhì)量分數(shù)Fig.4 Mass fraction of CO2in the non-hysteretic case over the entire model

圖5 滯后效應模型整個地層中各相態(tài)CO2質(zhì)量分數(shù)Fig.5 Mass fraction of CO2in the hysteretic case over the entire model

模擬到50年時，無滯后方案中氣相形態(tài)逐漸減少到70%，滯后方案中氣相保持約45%的比例。產(chǎn)生這種差異的原因是無滯后方案中，殘余氣相飽和度為0.177，即含水層介質(zhì)都默認捕獲此數(shù)值殘余氣相飽和度，只有氣相飽和度高于此殘余氣相飽和度時才運移。而在滯后方案中，是可能存在的最大氣相飽和度，即每個網(wǎng)格元素中根據(jù)經(jīng)歷的過程不同，轉(zhuǎn)換點飽和度不同，殘余氣相飽和度不同。因此滯后效應方案中捕獲的殘余氣相總是小于等于無滯后效應方案中的殘余氣相，因此表現(xiàn)出在模擬時間到50年時，無滯后效應方案中封存的氣相CO2比例大于滯后效應中的氣相比例。同樣，封存的CO2在含水層中的分布也不同，圖6為在CO2注入階段，即在疏干過程中氣相CO2飽和度的分布。在疏干過程中，兩種方案的相對滲透率特征曲線均沿主線變化，不涉及滯后現(xiàn)象，但因為殘余氣相飽和度的不同，相應特征曲線則不同。以注入點附近區(qū)域為例，無滯后效應的氣相krg總小于滯后效應方案的krg，因此在暈圖上表現(xiàn)出無滯后效應方案中氣相CO2運移緩慢(圖7)。

第4年末CO2停止注入后，在暈的尾端發(fā)生吸濕過程。如圖7，無滯后方案捕獲更多的殘余氣相，故CO2暈分布集中在運移過程中，在滯后方案中部分氣相CO2則快速運移到蓋層下，由于蓋層的封閉性好，氣相CO2暈逐漸向周圍擴散。

圖8表示數(shù)值計算中相對滲透率特征曲線。在無滯后方案中，特征曲線沿著疏干過程逆向變化，即相互驅(qū)替過程只有一條曲線。而在滯后方案中，相對滲透率和毛細壓力不僅與當前區(qū)域的飽和度有關(guān)，還和經(jīng)歷的過程和歷史飽和度有關(guān)。例如，在注入點網(wǎng)格區(qū)域，疏干過程中Sg增加，相應的krg增加，krl減少到0。在隨后停止注入的吸濕過程中，Sg逐漸減少，krg減少，krl逐漸增大。但在同一Sg時，兩個過程的相對滲透率不同，產(chǎn)生滯后的現(xiàn)象。同樣，在中間某網(wǎng)格處，相對滲透率變化趨勢相同，只是轉(zhuǎn)換點不同，進而相對滲透率變化曲線不同。

圖6 CO2注入過程中無滯后效應方案(上)和滯后效應方案(下)中氣相飽和度分布圖Fig.6 CO2plume evolution during the injection period in the non-hysteretic case(top)and the hysteretic case(bottom)

圖7 停止注入后無滯后效應方案(上)和滯后效應(下)方案中氣相飽和度隨時間變化分布圖Fig.7 CO2plume evolution during the post-injection period in the non-hysteretic case(top)and the hysteretic case(bottom)

圖8 模型中注入點網(wǎng)格相對滲透率變化曲線Fig.8 Relative permeability path for the injection grid

圖9～10對應兩種方案中毛細壓力特征曲線，選取模型中注入點網(wǎng)格和中間網(wǎng)格兩個不同位置得到其毛細特征曲線。滯后方案為一次完整疏干和吸濕過程，在中間網(wǎng)格處，氣相飽和度從0逐漸增大到轉(zhuǎn)換點飽和度，然后沿掃描曲線減少到相應殘余氣相飽和度，并未沿疏干過程主線變化，導致在同一氣相飽和度下毛細壓力不同。無滯后模型中，毛細壓力均沿主線變化。在實際情況中，僅用實驗測試得到的殘余氣相飽和度賦值給含水層來進行數(shù)值模擬，在封存能力、安全評估等方面會造成誤差。

圖9 滯后效應模型中不同位置毛細壓力變化曲線Fig.9 Capillary pressure path for grids in the hysteretic

圖10 無滯后效應模型不同位置毛細壓力變化曲線Fig.10 Capillary pressure path for grids in the non-hysteretic case

3 結(jié)論

在CO2咸水層封存的數(shù)值模擬計算中，相對滲透率和毛細壓力計算函數(shù)是重要計算參數(shù)之一。通過設(shè)計對應無滯后效應和考慮相對滲透率、毛細壓力滯后的兩種方案，運用數(shù)值模擬軟件TOUGH+CO2進行模擬，結(jié)果表明不同相對滲透率和毛細壓力計算函數(shù)的選擇表現(xiàn)出較大差異。滯后效應方案能更詳細描述每個網(wǎng)格中的相對滲透率和毛細壓力變化情況。在蓋層封閉性較好時，滯后效應對CO2封存總量沒有影響，但相態(tài)比例不同。

不同相對滲透率和毛細壓力計算函數(shù)的選擇表現(xiàn)出CO2暈運移分布不同，對實際封存方案設(shè)計，封存能力評估以及環(huán)境風險評價等都會帶來不同影響，因此在CO2咸水層封存數(shù)值模擬中，尤其是停止注入之后應引入與滯后現(xiàn)象相關(guān)的模擬計算。

［1］郝艷軍，楊頂輝.二氧化碳地質(zhì)封存問題和地震監(jiān)測研究進展［J］.地球物理學進展，2012，27(6):2369 － 2383.［HAO Y J，YANG D H.Research progress of carbon dioxide capture and geological sequestration problem and seismic monitoring research［J］.Progress in Geophys，2012，27(6):2369－2383.(in Chinese)］

［2］張煒，呂鵬.二氧化碳地質(zhì)封存中對流混合過程的研究進展［J］.水文地質(zhì)工程地質(zhì)，2013，40(2):101－107.［ZHANG W，LV P.Density-driven convection in carbon dioxide geological storage:a review ［J］.Hydrogeology ＆ Engineering Geology，2013，40(2):101－107.(in Chinese)］

［3］Juanes R，Spiteri E J，Orr F M，et al.Impact of relative permeability hysteresis on geological CO2storage［J］.Water Resource Research，2006，42(12):1－13.

［4］Doughty C.Investigation of CO2plume behavior for a large-scale pilot test of geologic carbon storage in a saline formation ［J］.Transport in porous media，2010，82(1):49 －76.

［5］許雅琴，張可霓，王洋.基于數(shù)值模擬探討提高咸水層CO2封存注入率的途徑［J］.巖土力學，2012，33(12):3825－3832.［XU Y Q，ZHANG K N，WANG Y.Numerical investigation for enhancing injectivity of CO2storage in saline aquifers［J］.Rock and Soil Mechanics，2012，33(12):3825 －3832.(in Chinese)］

［6］張曉宇，成建梅，劉軍，等.CO2地質(zhì)處置研究進展［J］.水文地質(zhì)工程地質(zhì)，2006，33(4):85－89.［ZHANG X Y，CHENG J M，LIU J，et al.An overview of underground sequestration of carbon dioxide［J］.Hydrogeology ＆ Engineering Geology，2006，33(4):85 －89.(in Chinese)］

［7］周琦，朱學謙，劉傳喜.毛細管力滯后曲線和相滲透率滯后曲線及處理［J］.油氣采收率技術(shù)，1999，6(3):47－50.［ZHOU Q，ZHU X Q，LIU C X.Curves of hysteretic phenomenon in relative permeability and capillary pressure［J］.Oil＆ Gas Recovery Technology，1999，6(3):47－50.(in Chinese)］

［8］Finsterle S.Multiphase inverse modeling:Review and iTOUGH2 applications ［J］.Vadose Zone Journal，2004，3(3):747－762.

［9］Zhang Y Q，Oldenburg C M，F(xiàn)insterle S，et al.Systemlevel modeling for economic evaluation of geological CO2storage in gas reservoirs［J］.Energy Conversion and Management，2007，48(6):1827 －1833.

［10］Zhang K，Doughty C，Wu Y S，et al.Efficient parallel simulation of CO2geologic sequestration in saline aquifers［C］//Lawrence Berkeley National Laboratory.SPE Reservoir Simulation Symposium.Berkeley:Society of Petroleum Engineers，2007.

［11］Yamamotoa H，Zhang K，Karasakib K，et al.Largescale numerical simulation of CO2geologic storage and its impact on regional groundwater flow:A hypothetical case study at Tokyo Bay［J］.Japan Energy Procedia，2009，1(1):1871－1878.

［12］Pruess K，Oldenburg C，Moridis G.TOUGH2 User’s Guide［R］. LBNL-43134.Berkeley:Lawrence Berkeley National Laboratory，1999.

［13］Zhang K，Moridis G，Pruess K.TOUGH+CO2:A multiphase fluid-flow simulator for CO2geologic sequestration in saline aquifers［J］.Computers ＆Geosciences，2011，37(6):714 －723.

［14］Finsterle S.iTOUGH2 user’s guide［R］.LBNL－40040.Berkeley:Lawrence Berkeley National Laboratory，1999.

［15］Doughty C. Modeling geologic storage of carbon dioxide Comparison of non-hysteretic and hysteretic characteristic curves［J］.Energy Conversion and Management，2007，48(6):1768 －1781.

［16］王大純，張人權(quán)，史毅虹，等.水文地質(zhì)學基礎(chǔ)［M］.北京:地質(zhì)出版社，1995.

［17］薛禹群.地下水動力學［M］.北京:地質(zhì)出版社，1997.

［18］Van Genuchten M Th，A closed-form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Science Society of America Journal，1980，44(5):892－898.