鄂爾多斯盆地深部咸水層二氧化碳地質(zhì)儲(chǔ)存熱-水動(dòng)力-力學(xué)(THM)耦合過程數(shù)值模擬

雷宏武，李佳琦，許天福，王福剛

吉林大學(xué)地下水資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)春 130021

?

鄂爾多斯盆地深部咸水層二氧化碳地質(zhì)儲(chǔ)存熱-水動(dòng)力-力學(xué)(THM)耦合過程數(shù)值模擬

雷宏武，李佳琦，許天福，王福剛

吉林大學(xué)地下水資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)春 130021

注入CO2到深部咸水層(CO2地質(zhì)儲(chǔ)存)被認(rèn)為是一種直接有效地減少CO2向大氣排放的途徑。CO2地質(zhì)儲(chǔ)存涉及到熱、水動(dòng)力和力學(xué)耦合過程，該耦合過程是預(yù)測(cè)CO2在儲(chǔ)層中的遷移轉(zhuǎn)化、評(píng)價(jià)儲(chǔ)層儲(chǔ)存能力和分析潛在風(fēng)險(xiǎn)的關(guān)鍵?；赥erzaghi固結(jié)理論，在熱-水動(dòng)力(TH)耦合軟件TOUGH2框架中加入了力學(xué)模塊，形成了新的熱-水動(dòng)力-力學(xué)(THM)模擬器。結(jié)合鄂爾多斯盆地CO2捕獲和儲(chǔ)存(CCS)示范工程場(chǎng)地的地質(zhì)、水文地質(zhì)條件，采用新的THM模擬器數(shù)值分析了CO2注入后地層中的溫度、壓力、CO2飽和度、位移和有效應(yīng)力的時(shí)空變化特征。結(jié)果顯示：在井口保持8 MPa和35 ℃情況下，能夠?qū)崿F(xiàn)10萬(wàn) t/a的CO2注入量；壓力上升的范圍遠(yuǎn)遠(yuǎn)大于CO2運(yùn)移和溫度降低的范圍，注入20 a后，其最大距離分別達(dá)到接近邊界10 km、620 m和100 m；位移和應(yīng)力變化主要與壓力變化相關(guān)，注入引起最大抬升為0.14 m，在注入井附近位置儲(chǔ)層中有效應(yīng)力變化水平方向要大于垂直方向，而在遠(yuǎn)井位置相反；注入引起井附近有效應(yīng)力明顯減小，從而導(dǎo)致了孔隙度和滲透率的增大，增強(qiáng)了CO2注入能力。

CO2地質(zhì)儲(chǔ)存；熱-水動(dòng)力-力學(xué)耦合過程；數(shù)值模擬；鄂爾多斯盆地

0 前言

近年來(lái)全球氣候問題逐漸加劇，政府間氣候變化專門委員會(huì)(IPCC)的系列研究指出溫室氣體(尤其是CO2氣體)的大量排放是主要影響因素之一[1-2]，該結(jié)論已經(jīng)得到了世界范圍科學(xué)研究者的認(rèn)可，因此CO2減排已迫在眉睫。深部咸水層由于其分布廣、經(jīng)濟(jì)價(jià)值小和儲(chǔ)存潛力巨大而被公認(rèn)為是進(jìn)行規(guī)?；疌O2地質(zhì)儲(chǔ)存最有潛力的場(chǎng)所[1,3]。

注入大量CO2到咸水層中涉及到復(fù)雜的熱-水動(dòng)力-力學(xué)耦合過程，而該過程的刻畫是預(yù)測(cè)CO2的遷移轉(zhuǎn)化規(guī)律、評(píng)價(jià)儲(chǔ)存能力和分析其泄漏風(fēng)險(xiǎn)的關(guān)鍵。熱-水動(dòng)力-力學(xué)過程之間的互相耦合是異常復(fù)雜的，實(shí)驗(yàn)只能研究其中的部分耦合過程，而數(shù)值模擬能夠充分考慮它們之間的互相耦合過程。在CO2地質(zhì)儲(chǔ)存的耦合數(shù)值模擬中，研究最多的是水動(dòng)力和力學(xué)之間的耦合，對(duì)應(yīng)地使用最多的模擬工具是TOUGH2-FLAC3D[4]。該軟件由Rutqvist于2002年通過耦合熱-水動(dòng)力(TH)軟件TOUGH2和力學(xué)分析軟件FLAC3D形成，并成功應(yīng)用于深部咸水層注入CO2后CO2運(yùn)移和壓力上升[4]、蓋層力學(xué)變化[5]、最大CO2持續(xù)注入壓力[6]和Salah場(chǎng)地CO2注入后的地表位移分析[7]。最近TOUGHREACT和FLAC3D被聯(lián)合在一起用于研究CO2注入后的熱-水動(dòng)力-力學(xué)和化學(xué)耦合過程[8-9]。Winterfeld等[10]在TOUGH2的基礎(chǔ)上開發(fā)了并行的熱-水動(dòng)力和力學(xué)模擬器。Hu等[11]基于平均應(yīng)力法耦合TOUGH2建立了熱-水動(dòng)力-力學(xué)模擬器TOUGH2-EGS，并應(yīng)用到干熱巖中。

雖然TOUGH2-FLAC3D能夠很好地模擬CO2地質(zhì)儲(chǔ)存過程的熱-水動(dòng)力-力學(xué)(THM)耦合過程，然而兩個(gè)獨(dú)立的代碼使得其間的耦合是松散的，計(jì)算效率較低；同時(shí)FLAC3D代碼的封閉模式使得對(duì)力學(xué)模型刻畫的理解和改進(jìn)幾乎不可能。筆者首先基于Terzaghi固結(jié)理論和TOUGH2的熱-水動(dòng)力耦合模型，建立了熱-水動(dòng)力-力學(xué)耦合模型；然后在TOUGH2代碼的基礎(chǔ)上，開發(fā)了力學(xué)模塊形成了新的熱-水動(dòng)力-力學(xué)模擬器，并通過與解析解對(duì)比驗(yàn)證了其正確性；最后通過改進(jìn)的程序研究了鄂爾多斯盆地深部咸水層二氧化碳地質(zhì)儲(chǔ)存過程中的熱、水動(dòng)力和力學(xué)變化過程。

1 THM模擬器和驗(yàn)證

1.1 THM過程數(shù)學(xué)模型

地下多孔介質(zhì)中的流體(包括地層水和CO2)流動(dòng)一般符合達(dá)西定律，熱流包括對(duì)流和傳導(dǎo)兩部分?；谫|(zhì)量和能量守恒，可建立多相流動(dòng)和熱對(duì)流傳導(dǎo)方程[12]：

流動(dòng)過程為

熱對(duì)流傳導(dǎo)為

而在力學(xué)模型的建立過程中，由于地表是位移自由面，加之地層具側(cè)向無(wú)限延伸性，可以假設(shè)CO2注入過程中引起的垂直總應(yīng)力不變[6]，地層變形側(cè)向受限，只在垂直方向上下變形。由此得到有效應(yīng)力變化和垂向應(yīng)變方程。

有效應(yīng)力變化為

平均有效應(yīng)力變化為

垂向應(yīng)變?yōu)?/p>

TOUGH2系列[12]是國(guó)際上先進(jìn)的TH耦合過程數(shù)值模擬程序，其能夠刻畫多個(gè)過程和可信的模擬結(jié)果[13]，因此被廣泛用于CO2地質(zhì)儲(chǔ)存的研究中[5,14-16]。ECO2N[17]是TOUGH2的一個(gè)模塊，能夠刻畫H2O、CO2和NaCl系統(tǒng)中的多相流動(dòng)和熱對(duì)流傳導(dǎo)過程。為了充分利用它的優(yōu)勢(shì)，本文的力學(xué)模型直接被嵌入到該程序中。繼承TOUGH2的時(shí)空離散特征，本文改進(jìn)的THM模擬器(THM-Terzaghi)采用隱式時(shí)間差分和空間任意形狀的積分有限差。

1.2 THM過程耦合方法

通常情況下，THM過程耦合方法可以分為兩類：全耦合和部分耦合。采用全耦合方法可以得到準(zhǔn)確的結(jié)果，但將花費(fèi)大量的計(jì)算機(jī)時(shí)間，而部分耦合雖然結(jié)果精度上有所降低，但是算法實(shí)現(xiàn)相對(duì)容易。TOUGH2中TH的耦合采用的是全耦合法，而添加的力學(xué)采用的是部分耦合法，即在計(jì)算完溫度和壓力后直接利用力學(xué)模型計(jì)算有效應(yīng)力變化和應(yīng)變，并根據(jù)有效應(yīng)力的變化修改孔隙度和滲透率[4-5]：

式中：φr為殘余孔隙度；φ0為零應(yīng)力狀態(tài)下的孔隙度；k0為零應(yīng)力狀態(tài)下的滲透率；a和c為實(shí)驗(yàn)參數(shù)?？紫抖群蜐B透率的變化在Newton-Raphson迭代內(nèi)部過程中考慮，這樣可以盡快反饋到流動(dòng)過程中。

1.3 THM模擬器驗(yàn)證

為了驗(yàn)證開發(fā)的THM-Terzaghi程序的正確性，筆者通過與兩個(gè)存在解析解問題的對(duì)比進(jìn)行了驗(yàn)證。一維固結(jié)沉降問題(圖1，表1)的對(duì)比結(jié)果(圖2)顯示，THM-Terzaghi程序的計(jì)算結(jié)果和解析解差別不大，從而證明了THM-Terzaghi在水動(dòng)力-力學(xué)(HM)耦合方面具可靠性；一維溫度傳導(dǎo)問題(圖3，表2)的對(duì)比結(jié)果(圖4)顯示，THM-Terzaghi在溫度-力學(xué)(TM)耦合方面具可靠性。

圖1 一維固結(jié)沉降概念模型Fig.1 Conceptual model of 1D settlement induced by consolidation

圖2 一維固結(jié)沉降解析解和THM-Terzaghi數(shù)值解結(jié)果對(duì)比Fig.2 Comparison of analytical and numerical solutions for 1D settlement induced by consolidation

1.3.1 一維固結(jié)沉降解析解和數(shù)值解對(duì)比

對(duì)于上面的一維固結(jié)沉降問題，可以得到超孔隙水壓力分布和沉降的解析解[18]。

超孔隙水壓力為

地表位移為

式中：d為埋深；PL為荷載；H為厚度；TV為時(shí)間因數(shù)，TV=kEst/(μH2)，μ為水的黏度；Es為壓縮模量Es=Kv(1+ν)/[3(1-ν)]；n為自然數(shù)。

表1 一維固結(jié)沉降模型參數(shù)

Table 1 Model parameters for 1D settlement induced by consolidation

參數(shù)取值厚度/m50滲透率/m21.0×10-14黏度/(kg/(m·s))8.9×10-3體積模量/Pa8×107泊松比0.20Biot系數(shù)1.0荷載/(105Pa)3.0

表2 一維熱傳導(dǎo)引起沉降模型參數(shù)

Table 2 Model parameters for 1D settlement induced by thermal conduction

參數(shù)取值厚度/m50滲透率/m20.0密度/(kg·m-3)2550孔隙度0.01熱傳導(dǎo)系數(shù)/(W/(m·℃))2.34熱膨脹系數(shù)/K-11.5×10-6比熱容/(J/(kg·℃))690泊松比0.20初始溫度/℃60邊界溫度/℃10

Tb為上部邊界溫度；Ti為初始溫度。圖3 一維熱傳導(dǎo)引起的沉降概念模型Fig.3 Conceptual model of 1D settlement induced by thermal conduction

圖4 一維溫度傳導(dǎo)引起沉降解析解和THM-Terzaghi數(shù)值解結(jié)果對(duì)比Fig.4 Comparison of analytical and numerical solutions for 1D settlement induced by thermal conduction

1.3.2 一維熱傳導(dǎo)引起的沉降解析解和數(shù)值解對(duì)比

對(duì)于上面的一維熱傳導(dǎo)引起的沉降問題，可以得到溫度分布和沉降的解析解[19]。

溫度為

地表位移為

2 模型建立

2.1 鄂爾多斯CCS示范工程簡(jiǎn)介

鄂爾多斯盆地是我國(guó)陸上第二大沉積盆地，盆地面積約為24萬(wàn)km2，形成于晚三疊世。盆地內(nèi)含有豐富的石油、煤炭、天然氣和煤層氣等資源，是我國(guó)重要的能源基地。區(qū)內(nèi)存在或即將建成多個(gè)煤炭轉(zhuǎn)化等能源基地，包括煤制油、煤制甲醇、煤制烯烴和煤制天然氣等項(xiàng)目，這些煤化工項(xiàng)目建成后將形成大量的CO2集中排放源。目前已投產(chǎn)的煤化工項(xiàng)目每年排放CO2約870 萬(wàn)t；據(jù)初步估計(jì)，在建項(xiàng)目全部投產(chǎn)后將每年排放CO2近5 400 萬(wàn)t[20]。如此多的CO2排放將是制約當(dāng)?shù)孛夯ろ?xiàng)目的瓶頸之一。

為應(yīng)對(duì)全球氣候變化，實(shí)現(xiàn)CO2減排，神華集團(tuán)2002年以來(lái)陸續(xù)開展了一系列CO2捕獲和儲(chǔ)存(CCS)的研究工作，并于2010年開始正式實(shí)施中國(guó)第一個(gè)CCS示范工程建設(shè)，2011年1月試注成功，2013年4月累計(jì)注入10 萬(wàn)t。

注入場(chǎng)地位于伊金霍洛旗，離神華煤制油公司直線距離約為12km。注入井深約2 500m，依次穿過第四系，白堊系志丹群，侏羅系安定組、直羅組和延安組，三疊系延長(zhǎng)組、紙坊組、和尚溝組和劉家溝組，以及二疊系石千峰組、石盒子組、山西組，石炭系太原組、本溪組和奧陶系馬家溝組。整套地層中存在多套適于CO2地質(zhì)儲(chǔ)存的儲(chǔ)蓋層組合，其中石千峰和劉家溝儲(chǔ)蓋層組合被認(rèn)為是儲(chǔ)層條件最好和最具有潛力的層位[21]。

注入井共射孔18層，累計(jì)厚度為88m，其中石千峰組為52m，占60%。在離注入井50m和80m處各有一口監(jiān)測(cè)井，用于監(jiān)測(cè)壓力、溫度和CO2氣體飽和度等。采用多層單井混合方式注入，分層監(jiān)測(cè)。

2.2 概念模型

根據(jù)場(chǎng)地巖石的巖石礦物分析，儲(chǔ)層以含長(zhǎng)石類的砂巖為主，砂巖中石英和長(zhǎng)石體積分?jǐn)?shù)最大分別可占到50%和40%。礦物捕集過程中反應(yīng)緩慢，在百年的時(shí)間尺度內(nèi)，對(duì)CO2的礦物捕獲幾乎沒有貢獻(xiàn)。因此，在本文模型建立中忽略化學(xué)反應(yīng)過程，僅僅考慮了溫度、壓力和力學(xué)之間的耦合過程。

石千峰和劉家溝儲(chǔ)蓋層組合(圖5)是本區(qū)CO2地質(zhì)儲(chǔ)存最有潛力的層位，根據(jù)測(cè)井?dāng)?shù)據(jù)，初步估計(jì)其可容納注入CO2的一半左右。為了簡(jiǎn)化模型，筆者以該套組合為研究對(duì)象，分析CO2注入到地層后的溫度場(chǎng)、壓力場(chǎng)、位移和有效應(yīng)力，以及CO2暈運(yùn)移的時(shí)空演化過程。

圖5 石千峰和劉家溝儲(chǔ)蓋層組合地層特征描述Fig.5 Description of Shiqianfeng-Liujiagou storage-caprock formations

石千峰和劉家溝儲(chǔ)蓋層可用的儲(chǔ)層共有8層，累計(jì)厚度為51.8m，單層的最大厚度為9.0m，最小為4.4m。砂巖儲(chǔ)層的滲透率在10-15m2數(shù)量級(jí)上，而泥巖的在10-17m2數(shù)量級(jí)上；砂巖的孔隙度約為0.10，泥巖的約為0.03(圖5)。根據(jù)這些特征可知，石千峰儲(chǔ)層是低孔低滲的薄砂層儲(chǔ)層，劉家溝的超低滲泥巖是理想蓋層。

鄂爾多斯深部地層傾斜度為1°～3°，接近水平。對(duì)于水平分布的地層，THM模型可以采用以注入井為中心的徑向模型(圖6)來(lái)大大縮減模型中網(wǎng)格的數(shù)目。模型垂向上為考慮蓋層的影響，上下范圍分別擴(kuò)展到泥巖中，最終埋深為1 576～2 160m，由于泥巖的阻隔作用設(shè)定其為隔水邊界；水平方向上取10km位置給一類定壓力邊界條件。在CO2注入過程中，溫度的對(duì)流傳導(dǎo)較慢，影響范圍遠(yuǎn)遠(yuǎn)小于壓力的影響范圍，因此，上下邊界采用隔熱，側(cè)向邊界溫度恒定(實(shí)際上采用隔熱并不影響計(jì)算結(jié)果)。初始溫度根據(jù)Ti=10.5+0.031 9d(d為埋深)確定，初始?jí)毫Ψ撵o水壓力，初始垂向應(yīng)力假設(shè)僅僅是由于上覆地層的自重引起，初始水平應(yīng)力滿足側(cè)向受限的靜止應(yīng)力狀態(tài)。

為了最大限度地利用儲(chǔ)層的空間，在注入前CO2被壓縮到超臨界狀態(tài)(為了簡(jiǎn)化，本文稱為氣相) (P>7.382 MPa，T>31.04 ℃)。井筒中的流體流動(dòng)速度很快，遠(yuǎn)遠(yuǎn)大于其在儲(chǔ)層中的流動(dòng)；而相對(duì)孔隙介質(zhì)，井筒要光滑得多，因此，可以假設(shè)流體在其中的流動(dòng)能量損失忽略。那么，井筒中的壓力僅受重力影響，其滿足方程Pn+1=Pn+ρngΔz(Pn為井筒中第n個(gè)網(wǎng)格中的壓力；ρn為第n個(gè)網(wǎng)格的CO2密度；Δz為第n個(gè)網(wǎng)格的垂向長(zhǎng)度。Pn+1是第n個(gè)網(wǎng)格壓力和密度的函數(shù))。由于CO2的可壓縮性，其密度在井筒中并不恒定，是溫度和壓力的函數(shù)，需要采用“由上至下”的計(jì)算方法得到與儲(chǔ)層相連位置處的注入壓力。根據(jù)實(shí)際注入情況，井口注入溫度和壓力分別為35.0 ℃和8.0 MPa。

圖6 CO2注入二維徑向概念模型Fig.6 2D radial conceptual model for CO2 injection into aquifers

2.3 模型參數(shù)

地層基本參數(shù)見表3，相對(duì)滲透率和毛細(xì)壓力均采用Van Genuchten模型，泥巖的毛細(xì)進(jìn)入壓力比砂巖大兩個(gè)數(shù)量級(jí)。

在考慮力學(xué)的耦合效應(yīng)時(shí)，最為關(guān)鍵的參數(shù)是耦合方程中的實(shí)驗(yàn)參數(shù)a和c，這兩項(xiàng)根據(jù)前人[4-5]總結(jié)分別取5×10-8Pa-1和22.2。為了使平均有效應(yīng)力、孔隙度，以及滲透率在初始狀態(tài)保持一致，需要得到每一個(gè)位置的殘余和零應(yīng)力狀態(tài)的孔隙度，以及零應(yīng)力狀態(tài)下的滲透率。假設(shè)孔隙度變化為實(shí)際孔隙度的10%，即φ0-φr=0.1φ。由于砂巖

表3 二維徑向模型地層基本物理參數(shù)

儲(chǔ)層均為薄砂層，有效應(yīng)力在每一層中變化不大，因此，可以簡(jiǎn)化為計(jì)算每一儲(chǔ)層的殘余和零應(yīng)力狀態(tài)的孔隙度，以及零應(yīng)力狀態(tài)下的滲透率，計(jì)算結(jié)果見表4，其中忽略泥巖孔滲變化。

表4 孔隙度和滲透率變化模型參數(shù)

Table 4 Parameters for stress-depended porosity and permeability

儲(chǔ)層編號(hào)厚度/m孔隙度零應(yīng)力狀態(tài)殘余實(shí)際滲透率/(10-15m2)零應(yīng)力狀態(tài)實(shí)際R19.00.1050.0950.107.362.80R25.40.1260.1130.1215.405.70R35.60.1050.0950.104.521.60R47.60.0520.0470.050.280.10R55.40.1050.0950.105.181.80R64.40.1050.0950.106.942.40R77.20.1160.1040.1110.203.50R87.20.1370.1230.1319.606.60泥巖467.90.0300.0300.030.080.08

2.4 網(wǎng)格剖分

對(duì)于2D的徑向模型，網(wǎng)格由以注入井為中心的一系列圓環(huán)組成。單元在水平方向上共有101個(gè)，垂向上共有85個(gè)，合計(jì)8 585個(gè)。水平方向上網(wǎng)格大小由近井的0.2 m逐漸對(duì)數(shù)增長(zhǎng)到邊界的800.0 m。垂向上為了刻畫CO2的重力翻轉(zhuǎn)現(xiàn)象，儲(chǔ)層及附近進(jìn)行了網(wǎng)格加密，最小網(wǎng)格大小為2.4 m，最大為泥巖中的10.0 m(圖7)。

圖7 二維徑向網(wǎng)格剖分Fig.7 Grid generation for 2D radial model

3 模擬結(jié)果

3.1 溫度時(shí)空變化特征

相對(duì)冷的CO2注入到咸水層后將導(dǎo)致地層溫度的降低。砂巖儲(chǔ)層中熱以對(duì)流和傳導(dǎo)兩種形式進(jìn)行傳遞，而上下相鄰的泥巖中熱主要以傳導(dǎo)為主，故導(dǎo)致溫度主要沿著咸水層變化。注入1 a后，CO2注入引起的溫度變化在側(cè)向上的最大距離約為29 m，20 a后約為74 m(圖8)。在垂向上，注入20 a后部分層位的溫度變化已經(jīng)突破咸水層之間的泥巖隔層。與CO2暈的運(yùn)移和壓力影響距離相比，溫度的影響范圍要小很多。

3.2 壓力時(shí)空變化特征

對(duì)于側(cè)向?yàn)橐活惗▔毫吔绲哪Ｐ?，用于?chǔ)存注入的CO2的空間來(lái)自：有效應(yīng)力減小引起的多孔介質(zhì)的膨脹體積；地層水壓力增加導(dǎo)致的地層水壓縮體積；從一類邊界流出的地層水體積。在注入早期，壓力上升還未影響到邊界的情況下，CO2被儲(chǔ)存在前兩種空間里。對(duì)于固結(jié)成巖的多孔介質(zhì)，地層壓縮和流體壓縮能力(大約在10-10Pa-1數(shù)量級(jí)上)并不強(qiáng)，注入大量的CO2將引起大范圍的壓力上升。圖9顯示了兩種注入方案對(duì)R1、R3和R8儲(chǔ)層的壓力變化。根據(jù)井筒“由上至下”的計(jì)算方法，得到這3層的注入壓力分別為20.7、22.2和23.2 MPa。R1和R8的壓力差(23.2-20.7=2.5 MPa)稍微小于其初始狀態(tài)的壓力差(19.8-16.9=2.9 MPa)，主要是因CO2密度小于水引起的。壓力上升的側(cè)向范圍是注入1 a后2 km和10 a后達(dá)到7 km，很明顯早期壓力上升的范圍快于后期，主要是因?yàn)殡S著范圍的擴(kuò)大，注入的CO2被“平均”儲(chǔ)存到更大的空間中。

3.3 CO2飽和度時(shí)空變化特征

泥巖中較大的毛細(xì)壓力將阻止CO2穿透砂巖儲(chǔ)層上下的泥巖運(yùn)移，強(qiáng)制其沿儲(chǔ)層水平方向運(yùn)移。從圖10可以看到CO2的側(cè)向運(yùn)移最遠(yuǎn)距離在注入1 a和20 a后分別約為145 和616 m，遠(yuǎn)遠(yuǎn)小于壓力上升的范圍。CO2在第R1、R2和R8儲(chǔ)層中的運(yùn)移速度快，主要是因?yàn)槠渚哂休^大的滲透率。

3.4 垂向位移和有效應(yīng)力時(shí)空變化特征

垂向位移是由壓力和溫度的變化引起。注入冷的CO2到地層中會(huì)引起孔隙水壓力的增加和地層溫度的降低?？紫端畨毫Φ脑黾舆M(jìn)而引起有效應(yīng)力的減小，地層會(huì)向上隆起，而溫度的降低會(huì)導(dǎo)致地層的下沉，兩種相反的作用共同影響地層的垂向變形。從圖11a可以看到，地表最大隆起約為0.14 m，其位置并不在靠近井的地方，而是離井大概100 m，說(shuō)明在井附近溫度降低引起的地面下沉還是比較大的。從整個(gè)區(qū)域來(lái)看，整個(gè)注入期內(nèi)地面是隆起的，說(shuō)明地面變形還是壓力變化導(dǎo)致起主導(dǎo)作用，隆起范圍由注入1 a的2 km擴(kuò)展到20 a的接近10 km，這個(gè)范圍變化過程直接對(duì)應(yīng)于壓力上升的變化過程。

在水平方向上，由于側(cè)向受限，壓力和溫度變化均會(huì)引起水平有效應(yīng)力的變化；而在垂向上，地面的自由變形使得垂向有效應(yīng)力僅僅受壓力變化引起，而與溫度變化無(wú)關(guān)。圖11b顯示了井附近(x=0.3 m)和遠(yuǎn)井(x=2 000.0 m)位置注入20 a后的有效應(yīng)力變化?？梢钥吹?，近井位置垂向有效應(yīng)力在注入層位明顯降低，降低幅度為3～4 MPa，水平有效應(yīng)力在注入層位降低大于垂向有效應(yīng)力的變化，最大達(dá)到6.0 MPa；而在遠(yuǎn)井位置有效應(yīng)力整個(gè)垂直剖面上變化一致，垂直和水平有效應(yīng)力降低分別約為1.0 MPa和0.3 MPa。

3.5 注入速率

對(duì)于定壓力的注入方案，注入速率是場(chǎng)地儲(chǔ)存能力評(píng)價(jià)的一個(gè)重要指標(biāo)。采用注入井井口8.0 MPa和35.0 ℃的注入方案，CO2注入速率約為1.7 kg/s(等效5.4萬(wàn) t/a)，其中注入到儲(chǔ)層R8、R2、R1和R7的分別占到注入量的28%，21%，19%和14%，合計(jì)達(dá)到80%(圖12)。根據(jù)儲(chǔ)層厚度統(tǒng)計(jì)，可以近似認(rèn)為石千峰和劉家溝儲(chǔ)蓋層組合能夠容納一半注入的CO2，因此，10萬(wàn)t/a的CO2注入量采用此注入方案能夠?qū)崿F(xiàn)。

3.6 力學(xué)引起的滲透率變化

孔隙水壓力的增加和地層溫度的降低均引起有效應(yīng)力的減小，從而導(dǎo)致孔隙度和滲透率的增加。從圖13可以看到，整個(gè)區(qū)域的滲透率在CO2注入后均增加，其中在近井100 m內(nèi)增加最大，增加了10%～40%，這是因?yàn)樵诖朔秶鷫毫蜏囟茸兓^大。同時(shí)，也可以看到在最下層的儲(chǔ)層R8增加最大，主要?dú)w因于儲(chǔ)層R8的相對(duì)較高的初始孔隙度和滲透率使得進(jìn)入R8儲(chǔ)層的CO2較多，溫度和壓力變化也較大。

圖8 溫度變化時(shí)空分布特征Fig.8 Spatial and temporal distribution of change of temperature

圖9 壓力時(shí)空分布特征Fig.9 Spatial and temporal distribution of pressure

SG.CO2氣體飽和度。圖10 CO2飽和度時(shí)空變化特征Fig.10 Spatial and temporal distribution of CO2 saturation

圖11 垂向位移和有效應(yīng)力時(shí)空變化特征Fig.11 Spatial and temporal distribution of vertical displacement and effective stress

圖12 CO2注入速率時(shí)間變化Fig.12 Evolution of CO2 injection rate

圖13 力學(xué)引起的滲透率變化Fig.13 Change of permeability induce by mechanics

3.7 耦合過程對(duì)注入速率的影響

圖14顯示了考慮不同耦合效應(yīng)的注入速率。僅考慮水動(dòng)力過程(H)、耦合的熱-水動(dòng)力過程(TH)、耦合的水動(dòng)力-力學(xué)過程(HM)和耦合的熱-水動(dòng)力-力學(xué)過程(THM)的平均注入速率分別為1.66, 1.56, 1.78和1.76 kg/s，最大和最小之間相差約13%。對(duì)比H和TH與HM和THM，可以看到溫度效應(yīng)的影響：在不考慮力學(xué)效應(yīng)時(shí)，忽略溫度效應(yīng)后的差別為6%，而考慮力學(xué)效應(yīng)時(shí)，忽略溫度效應(yīng)的差別為3%。對(duì)比H和HM與TH和THM，可以看到力學(xué)效應(yīng)的影響：在不考慮溫度效應(yīng)時(shí)，忽略力學(xué)效應(yīng)后的差別為7%，而考慮溫度效應(yīng)時(shí)，忽略力學(xué)效應(yīng)的差別約為13%。由此可見，溫度對(duì)注入速率起負(fù)作用，而力學(xué)起正作用，力學(xué)正作用大于溫度負(fù)作用。深層的原因是溫度降低而導(dǎo)致流動(dòng)黏滯系數(shù)增加，增加了注入阻力，而有效應(yīng)力的減小增加了滲透率，減小了注入阻力。

圖14 耦合過程對(duì)注入速率的影響對(duì)比Fig.14 Comparison of total injection rate for the four models

4 結(jié)論

1) 基于Terzaghi固結(jié)沉降的理論和TOUGH2中的TH耦合模型，建立了THM耦合數(shù)學(xué)模型，并在TOUGH2的框架中加入了該力學(xué)過程，通過與解析解的對(duì)比驗(yàn)證了該程序的科學(xué)性。

2) 在井口注入壓力為8.0 MPa、溫度為35.0 ℃時(shí)，10 萬(wàn)t/a的注入量基本能夠得到保證，其中進(jìn)入到石千峰第R8、R2、R1和R7層的CO2占到總注入量的80%。

3) 通過THM數(shù)值模擬分析了CO2注入后儲(chǔ)層中溫度、壓力、CO2飽和度、位移和有效應(yīng)力的時(shí)空演化特征：壓力影響范圍最大，注入20 a基本影響到邊界；CO2運(yùn)移的距離其次，注入20 a最大為620 m；溫度影響范圍最小，注入20 a不到100 m；位移和應(yīng)力變化主要與壓力變化相關(guān)，注入引起最大抬升為0.14 m，在注入井附近位置儲(chǔ)層中有效應(yīng)力變化水平方向要大于垂直的，而在遠(yuǎn)井位置相反。另外，注入引起井附近有效應(yīng)力明顯減小，從而導(dǎo)致了孔隙度和滲透率的增大，增強(qiáng)了CO2注入能力。

4) 忽略力學(xué)過程將低估CO2進(jìn)入到儲(chǔ)層，而忽略溫度降低將高估其注入，力學(xué)的影響明顯大于溫度的影響。

[1] IPCC (Intergovernmental Panel on Climate Change). Carbon Dioxide Capture and Storage[R]. New York: Cambridge University Press, 2005.

[2] IPCC (Intergovernmental Panel on Climate Change). Climate Change 2007: The Physical Science Basis: Contribution of Working Group I to the Fourth Assessment Report of the IPCC[R]. New York: Cambridge University Press, 2007.

[3] Bachu S,Adams J J.Sequestration of CO2in Geological Media in Response to Climate Change: Capacity of Deep Saline Aquifers to Sequester CO2in Solution[J]. Energy Conversion and Management, 2003, 44(20): 3151-3175.

[4] Rutqvist J, Wu Y S, Tsang C F, et al. A Modeling Approach for Analysis of Coupled Multiphase Fluid Flow, Heat Transfer, and Deformation in Fractured Porous Rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2002, 39: 429-442.

[5] Rutqvist J, Tsang C F. A Study of Caprock Hydromechanical Changes Associated with CO2Injection into a Brine Formation[J]. Environmental Geology, 2002, 42: 296-305.

[6] Rutqvist J, Birkholzer J, Cappa F, et al. Estimating Maximum Sustainable Injection Pressure During Geological Sequestration of CO2Using Coupled Fluid Flow and Geomechanical Fault-Slip Analysis[J]. Energy Conversion and Management, 2007, 48:1798-1807.

[7] Rinaldi A P, Rutqvist J. Modeling of Deep Fracture Zone Opening and Transient Ground Surface Uplift at KB-502 CO2Injection Well, In Salah, Algeria[J]. International Jounal of Greenhouse Gas Control, 2013, 12:155-167.

[8] Taron J, Elsworth D, Min K B. Numerical Simulation of Thermal-Hydrologic-Mechanical-Chemical Processes in Deformable, Fractured Media[J]. International Journal of Rock Mechanics & Mining Sciences, 2009, 46: 842-854.

[9] 于子望, 張延軍, 張慶，等. TOUGHREACT搭接FLAC3D算法[J]. 吉林大學(xué)學(xué)報(bào)：地球科學(xué)版, 2013, 43(1): 199-206. Yu Ziwang, Zhang Yanjun, Zhang Qing, et al. Algorithm of TOUGHREACT Links to FLAC3D[J]. Journal of Jilin University: Earth Science Edition, 2013, 43(1): 199-206.

[10] Winterfeld P H, Wu Y S. Parallel Simulation of CO2Sequestration with Rock Deformation in Saline Aquifers[C]. [S.l.]: SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 2011.

[11] Hu L T, Winterfeld P H, Fakcharoenphol P, et al. A Novel Fully-Coupled Flow and Geomechanics Model in Enhanced Geothermal Reservoirs[J]. Journal of Petroleum Science and Engineering, 2013, 107:1-11.

[12] Pruess K, Oldenburg C, George M. TOUGH2 User’s Guide, Version 2.0[R]. Berkeley: Earth Science Division, Lawrence Berkeley National Laboratory, University of California, 1999.

[13] Pruess K, Garcia J, Kovscek T, et al. Code Intercomparison Builds Confidence in Numerical Simulation Models for Geological Disposal of CO2[J]. Energy, 2004, 29: 1431-1444.

[14] Pruess K, Garcia J. Multiphase Flow Dynamics During CO2Disposal into Saline Aquifers[J]. Environmental Geology, 2002, 42: 282-295.

[15] Xu T, Sonnenthal E, Spycher N, et al. TOUGHREACT ：A Simulation Program for Non-Isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geological Media: Application to Geothermal Injectivity and CO2Geological Sequestration[J]. Computers & Geosciences, 2006, 32: 145-165.

[16] 許天福, 金光榮, 岳高凡, 等. 地下多組分反應(yīng)溶質(zhì)運(yùn)移數(shù)值模擬：地質(zhì)資源和環(huán)境研究的新方法[J]. 吉林大學(xué)學(xué)報(bào)：地球科學(xué)版, 2012, 42(5): 1410-1425. Xu Tianfu, Jin Guangrong, Yue Gaofan, et al. Subsurface 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.

[17] Pruess K. ECO2N: A TOUGH2 Fluid Property Module for Mixtures of Water, NaCl, and CO2[R]. Berkeley: Earth Science Division, Lawrence Berkeley National Laboratory, University of California, 2004.

[18] 方云, 林彤, 譚松林. 土力學(xué)[M]. 武漢: 中國(guó)地質(zhì)大學(xué)出版社, 2003. Fang Yun, Lin Tong, Tan Songlin. Soil Mechanics[M]. Wuhan: China University of Geosciences Press,2003.

[19] Jaeger J C, Cook, N G W, Zimmerman R W.Fundamentals of Rock Mechanics[M].Malden: Blackwell Publishing, 2007.

[20] 任相坤, 崔永君, 步學(xué)朋，等. 鄂爾多斯盆地CO2地質(zhì)封存潛力分析[J]. 中國(guó)能源, 2010, 32(1): 29-32. Ren Xiangkun, Cui Yongjun, Bu Xuepeng, et al. The Potential Analysis of CO2Geological Storage in Ordos Basin[J]. Energy of China, 2010, 32(1): 29-32.

[21] 吳秀章, 崔永君. 神華10萬(wàn)t/a CO2鹽水層封存研究[J].石油學(xué)報(bào): 石油加工, 2010(增刊1): 236-239. Wu Xiuzhang, Cui Yongjun. Study on 0.1 Mt/a CO2Sequestration in Saline Formation[J]. Acta Petrolei Sinica: Petroleum Processing Section, 2010(Sup.1): 236-239.

歡 迎 投 稿

《吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)》是由教育部主管、吉林大學(xué)主辦的地學(xué)類學(xué)術(shù)期刊。本刊主要刊發(fā)地質(zhì)與資源、地質(zhì)工程與環(huán)境工程、地球探測(cè)與信息技術(shù)等學(xué)科的原創(chuàng)論文，為地球科學(xué)領(lǐng)域中的新思想、新觀點(diǎn)、新理論、新方法、新技術(shù)和新應(yīng)用提供交流平臺(tái)。

本刊為中文核心期刊(GCGJ)、中國(guó)科技核心期刊(ISTIC)、中國(guó)科學(xué)引文數(shù)據(jù)庫(kù)重要期刊(CSCD)、中國(guó)學(xué)術(shù)期刊評(píng)價(jià)(2013-2014)權(quán)威期刊(RCCSE)。

本刊為“2012、2014中國(guó)最具國(guó)際影響力學(xué)術(shù)期刊”,連續(xù)4屆獲得教育部高校優(yōu)秀科技期刊獎(jiǎng)和第5屆教育部高校精品科技期刊獎(jiǎng),被國(guó)內(nèi)外SCOPUS,GEOREF,CA等20余家數(shù)據(jù)庫(kù)和文摘刊物收錄和摘引。

網(wǎng)上投稿：http://xuebao.jlu.edu.cn/dxb

Numerical Simulation of Coupled Thermal-Hydrodynamic-Mechanical (THM) Processes for CO2Geological Sequestration in Deep Saline Aquifers at Ordos Basin, China

Lei Hongwu, Li Jiaqi, Xu Tianfu, Wang Fugang

KeyLaboratoryofGroundwaterResourcesandEnvironment,MinistryofEducation,JilinUniversity,Changchun130021,China

Injecting CO2into deep saline aquifers, referred to as CO2geological sequestration (CGS), is considered to be a promising method to reduce the emission of anthropogenic CO2to atmosphere. CGS involves coupled thermal-hydrodynamic-mechanical (THM) processes, which are important for predicting migration and transformation of injected CO2, evaluating the reservoir performance, and assessing the risk associated. Based on the Terzaghi consolidation theory, a coupled mechanical module is developed that is incorporated into the existing simulator TOUGH2, which is a well-established code for TH processes in subsurface flow systems. Based on the geological and hydrogeological conditions of the Ordos CCS Demonstration Project, the new THM simulator is used to numerically analyze the spatial and temporal distribution of the temperature, pressure, CO2saturation, vertical displacement, and effective stress. The results show that 105metric tons of CO2per year can be finished under the injection condition of wellhead pressure 8 MPa and 35 ℃. The lateral distance with pressure buildup is more than that of CO2plume and that of temperature decrease . They are 10 km, 620 m and 100 m, respectively after 20 years’ injection. The vertical displacement and change of the effective stress are mainly related to the change of the pressure. The maximum surface uplift is about 0.14 m. The effective stress change is more significant in the horizontal direction than that in the vertical direction near the injection well, while it’s in the opposite away from the injection well. Injection induces an obvious decrease in the effective stress but enhances the porosity and permeability, this, in turn increases the CO2injectivity ultimately.

CO2geological sequestration; coupled thermal-hydrodynamic-mechanical processes; numerical simulation; Ordos basin

10.13278/j.cnki.jjuese.201502204.

2014-05-06

中國(guó)地質(zhì)調(diào)查局工作項(xiàng)目(12120113006300)；國(guó)土資源部公益性行業(yè)科研專項(xiàng)項(xiàng)目(201211063-06)；吉林大學(xué)博士交叉學(xué)科研究項(xiàng)目(2012JC014)

雷宏武(1985--)，男，博士，主要從事多相流數(shù)值模擬程序的開發(fā)和應(yīng)用研究，E-mail:hongwulei2008@aliyun.com。

10.13278/j.cnki.jjuese.201502204

P641.69；TK529

A

雷宏武，李佳琦，許天福，等.鄂爾多斯盆地深部咸水層二氧化碳地質(zhì)儲(chǔ)存熱-水動(dòng)力-力學(xué)(THM)耦合過程數(shù)值模擬.吉林大學(xué)學(xué)報(bào):地球科學(xué)版,2015,45(2):552-563.

Lei Hongwu, Li Jiaqi, Xu Tianfu, et al. Numerical Simulation of Coupled Thermal-Hydrodynamic-Mechanical (THM) Processes for CO2Geological Sequestration in Deep Saline Aquifers at Ordos Basin, China.Journal of Jilin University:Earth Science Edition,2015,45(2):552-563.doi:10.13278/j.cnki.jjuese.201502204.