应用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]按空间数据模型可分为面模型、体模型及混合模型三大类型。基于面模型构模是指将多个面元组合起来，通过表面形成三维空间目标轮廓，其重点在于空间实体的表面，如地形表面、岩层层面以及地下构筑物的轮廓；该方法的优点是数据存储量小，建模速度快，且便于显示和数据更新，缺点是不利于空间分析。体模型数据结构侧重于空间实体边界与内部的整体表示，如矿体、水体等，通过对体的描述实现三维空间目标表示；其优点是适于空间操作和分析，但运算量、存储空间占用大，构模速度慢。混合模型为面模型与体模型的综合。常见三维地质建模方法见表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.