應用TOUGH模擬二氧化碳地質儲存過程的復雜地質體建模技術與實現

楊艷林，許天福，李佳琦，王福剛

吉林大學地下水資源與環境教育部重點實驗室，長春 130021

應用TOUGH模擬二氧化碳地質儲存過程的復雜地質體建模技術與實現

楊艷林，許天福，李佳琦，王福剛

吉林大學地下水資源與環境教育部重點實驗室，長春 130021

基于TOUGH2數值模擬器中存在處理復雜地質條件難的問題，提出了將三維地質模型與TOUGH2數值模型耦合的具體思路和方法，并在Windows平臺上編制了相應的程序，利用該程序可以直接將地質模型(如GMS和Petrel建立)與TOUGH2的數值模型進行有機融合，并進行了算法驗證。通過CO2地質儲存的實例，可以看出地層構造對CO2的空間運移起著控制作用，模擬結果顯示CO2會沿著背斜面地層進行運移擴散。因此通過此轉換程序，能夠增強TOUGH2模擬器處理地層起伏、斷層、褶被等常見地質構造，以及地質屬性隨空間變化的地質體，提高了模擬器的使用效率。

地質模型；TOUGH2；數值模擬；CO2地質儲存

0 引言

隨著工程地質學、數學地質學、圖形學、計算機科學以及其他學科的快速發展和應用，數值模擬技術在工程領域中的應用越來越廣泛。二氧化碳地質儲存作為有效的二氧化碳減排手段，成為應對全球氣候變化問題研究中的熱點之一。由于二氧化碳地質封存實際操作與監測難度大、費用高昂，因此數值模擬作為一種有效而經濟的技術方法，成為二氧化碳地質儲存研究全過程的重要工具。

數值模擬結果的可靠性在很大程度上取決于對巖體基本性質的認識和各種地質因素的合理簡化。在二氧化碳地質儲存數值模擬領域，應用最為廣泛和成功的軟件為TOUHG[1-4]家族軟件系列。TOUGH2模擬器中的ECO2N模塊在這方面表現最為出色，但是其卻難以處理具有復雜地質構造(如斷層、巖脈、褶皺，圖1)以及巖石物性特征(如孔隙度、滲透率)空間變化的地質體；而融合了沉積學、構造地質學、地震地層學和測井地質學的專門性的地質三維模擬模型[6]，能較合理地描述研究區的構造、沉積微相以及物性參數的空間分布。如果能夠將已有的三維地質建模軟件與TOUGH2模擬軟件進行技術上的耦合，彌補TOUGH2模擬器在復雜地質體建模方面的不足，將大大提升TOUGH2軟件對于實際場地二氧化碳地質儲存模擬結果的精準度和可靠性，對于工程決策和安全性評價具有十分重要的意義。

黑箭頭為CO2的運移方向；白箭頭為咸水運移方向[5]。圖1 CO2泄露過程示意圖Fig. 1 Underground CO2 leakage

1 三維地質建模技術概述

三維地學建模(3D geosciences modeling，3DGM)，是首先由加拿大工程地質學家Simon W. Houlding[7-8]提出的。它主要是利用計算機和科學可視化技術，直接在3D空間中以數字化的形式表達和再現地質體與地質環境，進而輔助工程設計、施工與決策。經過十多年的研究和發展，三維地質模型研究已經在石油勘探、礦山開采、工程地質等領域廣泛應用，并研發了不少成熟的商業軟件，如美國Reservoir Characterization Research and Consulting公司的3D Earth Modeling軟件、中國石油大學開發的“RDMS”軟件等。

現有的地質建模方法[9-12]按空間數據模型可分為面模型、體模型及混合模型三大類型。基于面模型構模是指將多個面元組合起來，通過表面形成三維空間目標輪廓，其重點在于空間實體的表面，如地形表面、巖層層面以及地下構筑物的輪廓；該方法的優點是數據存儲量小，建模速度快，且便于顯示和數據更新，缺點是不利于空間分析。體模型數據結構側重于空間實體邊界與內部的整體表示，如礦體、水體等，通過對體的描述實現三維空間目標表示；其優點是適于空間操作和分析，但運算量、存儲空間占用大，構模速度慢?；旌夏Ｐ蜑槊婺Ｐ团c體模型的綜合。常見三維地質建模方法見表1。

2 TOUGH數值模擬器

TOUGH2[1]是美國勞倫斯伯克利國家實驗室開發的用于模擬一維、二維和三維孔隙或裂隙介質中多相流多組分非等溫流動過程的模擬器，它采用積分有限差分方法進行空間離散，運用全隱式的迭代技術對數值模型進行求解，被廣泛用于地熱、核廢物處置、二氧化碳地質儲存、環境污染評價和修復等方面的研究，是一款通用的地下流體數值模擬軟件。

對于任意區域Vn，它的質量(水、氣或化學組分)和能量方程可表述為

表1 常見三維地質建模模型方法分類

等式左邊：

對于β相，等式右邊可寫為

a. 某單元；b. 平面上相鄰兩單元。Fnm為截面積；Dn，Dm為距離；Anm為表面積。圖2 典型網格示意圖Fig. 2 Schematic representation of 2D grid

由于采用積分有限差進行空間離散，故只需關心單元的屬性(如體積、面積、單元所具有的屬性以及繪圖所用的單元點坐標，圖2a)以及與相鄰單元的連接信息(如連接的方向和與重力方向的夾角余弦值、接觸面積和單元點到接觸面的距離，圖2b)即可。基于此特點，常用的網格剖分有徑向網格、規則格網以及局部正交網格(perpendicular bisectors)(圖3)。由于徑向網格是對將三維問題二維化的簡化，不便于三維問題的研究，本次研究不對其進行地質模型耦合。

a, b. 規則格網模型以及相應的三維模型；c, d. 不規則格網模型以及相應的三維模型。圖3 網格模型示意圖Fig. 3 Sketched map of typical grid model

3 TOUGH2的復雜地質體三維地質建模技術與實現

在地下水模擬領域，將地質模型屬性直接轉化為水流數值模型已得到了廣泛的應用。如地下水數值模擬中，可視化界面GMS[13](groundwater model system)運用Solids模塊建立場地三維可視化地層模型，之后直接轉到地下水數值模型(MODFLOW[3]計算模型)中(Solids→MODFLOW)，可以逼真地刻畫地層的空間結構，能有效地處理復雜的地層構造特征(出露和尖滅)等問題；在油藏數值模擬中功能強大的Petrel軟件[14-16]可以為油藏數值模擬提供精細的三維地質模型，經網格粗化后轉為油藏數值模擬所需的油藏模型。

筆者擬對二氧化碳地質儲存TOUGH數值模擬器做一些嘗試性的工作，以增強其對復雜地質體的建模能力。

通過本文前面對地質模型與TOUGH模擬器的特點分析知，將地質模型屬性轉到數值模型中，主要涉及到2種網格模型的匹配操作，如對于以四面體為數據結構(表1)的地質模型和以矩形格網剖分的TOUGH數值模型，就是要對矩形格網單元與四面體單元進行匹配操作。一種簡單而實用的匹配方法是通過坐標來進行判定，如地質模型的單元是否包括TOUGH模擬器網格單元的坐標，若在其內，則繼承地質模型單元所具有的屬性;若在其上，則找出相鄰的單元，取相鄰單元的加權屬性。否則不能繼承地質模型單元的屬性。這一過程的大致流程為：首先讀入TOUGH2模型的單元坐標信息；然后將其與地質模型網格進行遍歷比較；最后確定單元屬性(其算法見圖4)。在這個過程中主要是判定點與地質模型單元的關系，依據前面對地質模型的分析，可將其總結為點與矩形網格、點與三棱柱(也包括似三棱柱)、點與四面體等關系的判定，其中都要涉及點與平面的關系，下面將其實現算法進行具體闡述。

3.1 點與平面的關系

在三維空間中，取截面上不共線的3點(圖5a)，記為(x1，y1，z1)，(x2，y2，z2),(x3，y3，z3)，則其確定的平面方程為

其中：

E1=z1(y3-y2)+z2(y1-y3)+z3(y2-y1) ；

E2=z1(x2-x3)+z2(x3-x1)+z3(x1-x2) ；

E3=y1(x3-x2)+y2(x1-x3)+y3(x2-x1) ；

E4=z1(y2x3-y3x2)+z2(y3x1-y1x3)+

z3(y1x2-y2x1) 。

其法向量為：n=(E1,E2,E3)(在計算時，保證E3≥0)。對于空間中的任意一點(x0,y0,z0)，令

F(x0,y0,z0)=E1x0+E2y0+E3z0+E4。

規定：若F(x0,y0,z0)>0，則點在面的“上方”；若F(x0,y0,z0)=0，則點在平面上；若F(x0,y0,z0)<0，則點在平面的“下方”。

3.2 點與網格單元的判定

對于柵格和八叉樹體元模型，在轉換的過程中需考慮點與矩形網格的關系，不考慮側面不共面的情況，因為這種網格單元地質建模較少，見圖5b。在判定的過程中主要有2個步驟：首先判定點在單元水平投影的矩形內；其次判定點在上下表面之間，即在面ABCD與平面A1B1C1D1之間。其中，由于A、B、C、D中4點可能不在同一平面內,故將其分解為2個平面來處理，即平面ABC與平面ACD(平面ABCD類似)，記為FABC=F(x0,y0,z0)(FABC為ABC的平面方程)與FACD=F(x0,y0,z0)。在運用判定時，需先確定點所在的平面是在平面ABC內還是在平面ACD內，后再計算對應的平面方程，并計算對應的FABCD值。若在平面ABC內，則計算FABCD=FABC，若FABCD=0，則點在其面上。同理，對面A1B1C1D1進行同樣的處理，并計算對應的值FA1B1C1D1，再計算Fgrid=FABCDFA1B1C1D1：若Fgrid<0，點在其間；Fgrid>0，則在其外。

圖4 算法耦合過程Fig. 4 Process of coupling algorithm

對于層面地質模型，一般用三棱柱作為單元(指上下為三角形，側面具有共面特征，在水平面上的投影為三角形，圖5c)。在判定點是否在三棱柱內，主要有2個過程：首先判定點是否在投影三角形內；然后再判定點是否在上下2個三角形內(與點和矩形網格單元的判定法類似)。

在通過鉆孔構建地質模型時，由于鉆孔開采時的諸多原因，導致鉆孔在垂直方向上會產生一定的偏移，使得上、下兩相鄰的表面TIN上對應點無法保證在同一鉛垂線上，因而導致所連接形成的三棱柱的側面不在一個平面內，即A、A1、B1、B，A、A1、C1、C或B、B1、C1、C每組中的4點不共面(圖5d)，這種三棱柱即為似三棱柱[4]。其判定較復雜，大致流程為：首先通過上下平面，確定點是否在三棱柱之間，即點O是否在平面ABC與平面A1B1C1之間；若在其間，然后判定點是否在3個側面內(ABA1B1、BCB1C1和CAC1A1)。由于各個面可能不共面，故將其分為2個面進行處理，類似于點與矩形網格中的第二步；再通過點與各個面的關系，確定點在側面內、外或上。

a. 點與平面；b. 點與格網單元；c,d. 點與三棱柱；e. 點與四面體。圖5 點與體的關系示意圖Fig. 5 Relation diagram of point and body

在真三維地質模型中，大部分地質單元是以四面體為數據結構的。在確定點與地質模型單元關系時，就是要確定出點所在的四面體單元(圖5e)。其完全可以通過點與相鄰面的關系得出點是在四面體內、外還是面上。以其中的一組面(平面ABD和平面CBD)為例：若FABD=0或FCBD=0，則點在面上；若FABDFCBD>0，則點在四面體外；若FABDFCBD<0，則需繼續判定點與其他組面的關系。在這一過程中，對于一個四面體單元最多需進行3次判定。

4 算例驗證

基于上述理論和方法，在Windows平臺上利用VC++和OpenGL開發工具，開發了TOUGH數值模擬軟件的轉換程序及可視化界面，解決了與其他地質建模軟件(如Petrel軟件，GMS的Solids模塊建模軟件等)耦合的技術問題，實現了直接將其地質體模型數據整合到TOUGH模型中。

利用某實際研究區的測試數據驗證了該技術的科學性。該研究區具有背斜的地質構造，利用GMS的Solids模塊建立了對應的地質模型，由于其是采用TIN法來進行構建的，故其基本單元是三棱柱單元，后通過中間程序判定點與三棱柱的關系，將其耦合到了TOUGH模型中。其中：圖6a為研究區建立地質模型的平面網格剖分情況；圖6b為研究區的地質模型和剖面切割圖；圖6c為TOUGH模型的網格剖分平面圖，其采用局部加密的矩形網格剖分方式；圖6d為其對應的概念模型，表層為蓋層，第2層為深部含水層，第3層為泥巖蓋層，第4層為CO2儲層，第5層為基巖(從上往下)；圖6e為利用開發的中間轉換程序，將生成地質模型耦合到了TOUGH2模型中的巖性圖以及剖面圖，與實際地質模型吻合較好；圖6f為CO2在80 a后的空間運移圖和空間等值面圖，可明顯看出，CO2沿著背斜面地層進行運移擴散。

a. 地質模型的平面網格圖；b. 3D地質模型圖與剖面圖；c. TOUGH平面網格模型；d. 概念模型；e. 地質模型轉到TOUGH數值模型的3D網格與剖面圖；f. 80 a后CO2飽和度的空間分布圖與等值面。圖6 耦合模型驗證Fig. 6 Verify the coupling model

5 結論與討論

1)將已有三維地質建模軟件的地質模型耦合到TOUGH2模擬器的數值模型中，以致所做的數值模型繼承了實際地質情況，使模擬結果更能反映實際；同時對于有各種地質構造(如地層起伏和斷層等)的地質體在數值模擬中的處理也得到了解決。

2)通過將其運用到一個進行CO2地質封存的背斜地質構造的實例中，可以看出CO2沿著背斜地層進行擴散運移，與實際相符合；同時從另一方面可知地質構造對CO2的空間運移起著控制作用，若不能對其進行正確的處理，即有可能得出不正確的結果，從而印證了該耦合算法的正確性和有效性。

3)在耦合算法部分，本文對其實現進行了詳細描述，但實際上，數值模型單元可能與多個地質模型單元相交，因此需要按照其相交的體積進行地質屬性加權平均，尤其是對于地質屬性變化較大的地區，這方面還需進行深入的研究。

[1] Pruess K,Oldenburg C,Moridis G. TOUGH2 User’s Guide, Version 2.0[R]. Berkeley: Lawrence Berkeley National Laboratory, 1999.

[2] Xu T, Pruess K. Modeling Multiphase Non-Isother-mal Fluid Flow and Reactive Geochemical Transport in Variably Saturated Fractured Rocks: 1：Methodology[J]. Americal Journal of Science, 2001, 301: 16-33.

[3] Xu T, 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[R]. Berkeley: Lawrence Berkeley National Laboratory, 2005.

[4] 許天福，金光榮，岳高凡，等. 地下多組分反應溶質運移數值模擬：地質資源和環境研究的新方法[J]. 吉林大學學報：地球科學版，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.

[5] Chang K, Bryant S L. The Effect of Faults on Dynamics of CO2Plumes4[J]. Energy Procedia, 2009,1(1):1839-1846.

[6] 程朋根，劉少華，王偉，等. 三維地質模型構建方法的研究及應用[J]. 吉林大學學報：地球科學版，2004，34(2)：309-313. Cheng Penggen, Liu Shaohua, Wang Wei, et al. Study and Application of a New 3D Geological Model Construction Method[J]. Journal of Jilin University：Earth Science Edition, 2004, 34(2): 309-313.

[7] Houlding S W. 3D Geoscience Modeling Computer Techniques for Geological Characterization4[M]. New York: Springer-Verlag, 1994.

[8] 鄔倫，劉瑜，張晶，等. 地理信息系統：原理、方法和應用[M]. 北京：科學出版社，2002. Wu Lun, Liu Yu, Zhang Jing, et al. Geographic Information System：Principle, Method and Application[M]. Beijing: Science Press, 2002.

[9] Li R. Data Structure and Application Issue in 3D Geographic Information System[J]. Gematics, 1994, 48(3): 209-244.

[10] 李清泉，李德仁. 三維空間數據模型集成的概念框架研究[J]. 測繪學報，1998，27(4):325-330. Li Qingquan, Li Deren. Research on the Conceptual Frame of the Integration of 3D Spatial Data Model[J]. Acta Geodaetica et Cartographica Sinica, 1998, 27(4):325-330.

[11] 李青元,朱小弟,曹代勇. 三維地質模型的數據模型研究[J]. 中國煤田地質,2000,12(3):59-63. Li Qingyuan, Zhu Xiaodi, Cao Daiyong. Data Model of 3D Geological Model[J]. Coal Geology of China, 2000, 12(3):59-63.

[12] Kreuseler M. Visualization of Geographically Related Multidimensional Data in Virtual 3D Scenes[J]. Computer and Geosciences, 2000, 26(1):101-108.

[13] Michael G M, Arlen W H. A Modular Three-Dimensional Finite-Difference Groundwater Flow Model[M]. Washington: United States Government Printing Office, 1988.

[14] 羅玉，翟中霞，孟慶友，等. 利用Petrel軟件進行精細地質建模研究[J]. 太原科技, 2008(10): 60-63. Luo Yu, Zhai Zhongxia, Meng Qingyou, et al. Research on Building Refined Geological Model Using the Software of Petrel[J]. Taiyuan Sci-Tech,2008(10):60-63.

[15] 向傳剛. PETREL軟件在油藏數值模擬研究中的應用[J]. 內蒙古石油化工,2009(15):8-11. Xiang Chuangang. Application of Petrel Software in the Study of Numerical Simulation of Oil Reservoir[J]. Inner Mongolia Petrochemical Industry, 2009(15): 8-11.

[16] Essca Group. Petrel 2005 Handbook of Software to Help[M]. New York：Schlumberger, 2004.

Complex Geological Body Modeling and Implementation of CO2Geological Storage Simulation Using TOUGH

Yang Yanlin，Xu Tianfu，Li Jiaqi，Wang Fugang

Key Laboratory of Groundwater Resources and Environment,Ministry of Education, Jilin University, Changchun 130021,China

It’s difficult to deal with complex geological conditions in TOUGH2 simulator, therefore, the method of a complex 3D geological model and coupling pattern numerical simulation were put forward, and the authors developed the corresponding program on the Windows platform. It can be used to the geological modeling created by GMS or Petrel software coupling TOUGH2 numerical model. Through the instance of carbon dioxide (CO2) geological storage, it can be seen that the stratigraphic structure plays a controlling role of CO2space migration, the simulation results show that CO2can spread along the back slope formation. Therefore, the coupled method,which provides a new approach and methodology to improve the efficiency of this simulator to deal with complex geological structure, such as ups and downs stratum, fault, fold and geological attribute of geological body chang according to the space.

geologic model; TOUGH2; numerical modeling; CO2geological storage

2013-11-26

中國地質調查局工作項目(12120113006300)

楊艷林(1984--)，男，博士研究生，主要從事多相流多組分數值模擬與程序開發方面的研究，E-mail:yang yanlinjida@gmail.com

王福剛(1975--)，男，副教授，主要從事水文地質與環境地質研究，E-mail:wangfugang@jlu.edu.cn。

10.13278/j.cnki.jjuese.201404209.

10.13278/j.cnki.jjuese.201404209

P66;P641

A

楊艷林，許天福，李佳琦，等. 應用TOUGH模擬二氧化碳地質儲存過程的復雜地質體建模技術與實現.吉林大學學報:地球科學版,2014,44(4):1307-1313.

Yang Yanlin，Xu Tianfu，Li Jiaqi, et al. Complex Geological Body Modeling and Implementation of CO2Geological Storage Simulation Using TOUGH.Journal of Jilin University:Earth Science Edition,2014,44(4):1307-1313.doi:10.13278/j.cnki.jjuese.201404209.