2. 中国科学院武汉岩土力学研究所, 武汉 430071
2. Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China
0 引言
GMS(groundwater modeling system)是美国杨伯翰大学环境模拟研究实验室和美国陆军排水工程试验站在大量数值分析代码(MODFLOW, MT3DMS, RT3D, SEAM3D, MPFPATH, MODAEM, SEEP2D, FEMWATER, WASH123D, UTCHEM)的基础上开发的一套用于地下水模拟前后处理的可视化界面[1-2]，由于界面友好、功能强大，在地下水领域得到了广泛的使用。TOUGH2是美国劳伦斯伯克利国家实验室开发的用于模拟一维、二维和三维孔隙和裂隙介质中多相多组分非等温流动过程的模拟器[3]，具有多个不同流体组合系统刻画的状态方程模块，被广泛用于地热、核废物处置、二氧化碳地质储存、环境污染评价和修复和非饱和带水文的研究等[4-5]，是一款通用的地下流体数值模拟软件。
GMS包含的MODFLOW(简写为GMS/MODFLOW)虽然只能用于地下浅层单相地下水的模拟，但是GMS具有强大的前后处理可视化界面[2, 6]。TOUGH2虽然数值分析功能强大并应用广泛，但是缺乏全面的前后处理功能，使得网格剖分、岩性属性赋值、初始边界条件设置和结果显示等模型建立和分析工作异常费时和困难。PetraSim是云端工程顾问公司(Thunderhead Engineering Consultants)专门为TOUGH家族软件开发的一套前后处理界面软件[7]；TOUGHer是犹他大学为TOUGHREACT开发的一个图形化界面软件[8]；pyTOUGH是采用python语言编写的用于TOUGH2三维模型建立的脚本库，由许多函数组成，直接调用可以自动产生TOUGH2所需文件[9]。这些软件都有很强的前后处理功能，但仍处于开发阶段。单独为TOUGH开发一个软件可能花费大量的时间，而通过联合TOUGH的数值计算能力和其他成熟软件的可视化功能是快速解决问题的捷径。前人在此方面进行了大量的研究：Audugane等[10]基于Petrel、Tecplot和ParaView编写了Petrel与TOUGH之间的互相转化程序，提供了一种快速建立TOUGH数值模型的方法；Borgia等[11]基于GMS产生的网格，联合TOUGH2进行了复杂地质环境下的三维地下水数值模拟，但没有直接考虑地层起伏情况；杨艳林等[12]将三维地质模型与TOUGH2数值模型结合，在Windows平台上开发复杂地质建模及TOUGH前后处理程序，使TOUGH2模拟器能处理地层起伏、断层及地质属性随空间变化等复杂的地质结构，Berry等[13]基于开源地理信息系统软件GRASS GIS，开发了TOUGHGIS前处理软件，增强了TOUGH2建模过程，最后模型计算结果使用TOUGH2 Viewer[14]进行显示；mView是一款商业前后处理软件[15]，兼容几乎所有的TOUGH2模块，能处理复杂地质体建模问题；Hu等[16]采用C++开发了前后处理程序，可基于平面泰森多边形建立复杂的三维模型。
本文在前人工作的基础上，主要进行了考虑地层起伏改进网格转化方法的研究，增加了GMS井文件和一类水头边界条件的转化；同时，还进行了整个GMS水流模型向TOUGH2/EOS1模型的转化，并通过两个实例对比了两个软件的数值模拟结果。
1 GMS/MODFLOW与TOUGH2原理对比及转化方法 1.1 GMS/MODFLOW和TOUGH2计算原理对比GMS/MODFLOW和TOUGH2最基本的计算原理都是质量和能量守恒，用于描述流体流动的均为达西定律，这样模型所需的参数在很大程度上是一样或类似的，从而为两个软件的对比提供了理论上的可行性；区别是，TOUGH2为了刻画多相流动行为，采用压力进行计算，相的物理化学属性参数均是随温度和压力变化的。具体为：GMS/MODFLOW仅将水头作为计算变量，其以单相流动且流动过程中密度不变为特征；而TOUGH2将压力、温度、饱和度、溶质质量分数作为计算变量，考虑多相流动和热对流扩散，流体属性(密度、黏度、热焓等)随着温度、压力和溶质浓度的变化而变化，不仅可以得到压力分布，而且可以得到相的饱和度和每个组分的质量分数。在原理上，两者的差异主要有：1) GMS/MODFLOW只有单一水组分；而TOUGH2可以考虑水、气(空气或CO2)和盐多种组分，因而可以模拟非包气带的流动，可以更加客观地模拟疏干问题。2) GMS/MODFLOW只适用于近地表的流动系统；而TOUGH2中的流体属性参数是温度和压力的函数，因此TOUGH2可以模拟高温高压体系中流体流动。3) GMS/MODFLOW不能模拟温度场；而TOUGH2能够模拟温度场，因此可适用于地热方面的分析和评价。4) GMS/MODFLOW内部组建的是线性方程组；而由于考虑属性变化，TOUGH2组建的是非线性方程组。针对承压含水层的模拟，GMS/MODFLOW的效率更高一些；而对于非线性非常强的问题，比如潜水含水层中的流动、承压水和非承压水之间的反复转化和非均质程度高的含水层等，TOUGH2的可靠性和效率更高一些。5) 在数值离散方面，GMS/MODFLOW采用的是规则六面体有限差分方法；而TOUGH2采用的是不规则单元的积分有限差，TOUGH2比GMS/MODFLOW更加通用和灵活。
GMS/MODFLOW中控制方程的一般形式[17]为
其中：Kxx，Kyy，Kzz分别为沿x，y，z方向的渗透系数，m/s；h为水头，m；W为单位时间单位体积源汇，s-1；Ss为孔隙介质的储水率，m-1；t为时间，s。
TOUGH2中控制方程的一般形式[3]为
流动过程为
热对流传导为
其中：M为质量或能量，kg或J；F为质量或能量通量，kg/m2或J/m2；n为边界外法向向量；Γ和Γn表示网格间的连接面和连接面面积，m2；V和Vn表示微元体和体积，m3；κ=w, i, g分别表示组分水、盐、气；κ+1表示热或能量；β=A, G分别表示气相和液相；uβ为β相的内能；φ为孔隙度；g为重力加速度，m/s2；μβ为液体的动力黏滞系数，Pa·S；Sβ为β相的饱和度；ρβ和ρR分别为β相和岩石的密度，kg/m3；X为质量分数；CR为岩石比热，J/(kg·K)；k和kr分别为渗透率和相对渗透率，m2；p为压力，Pa；T为温度，K；λ为导热系数，W/(m·K)；hβ为β相的焓值，J/kg；J为弥散项，kg/m2；q为单位时间源汇，kg/s。
1.2 GMS/MODFLOW数值模型文件转化为TOUGH2输入文件GMS/MODFLOW模型的建立有两种方法：概念模型法和网格法。概念模型法基于刻画模型的点、线、区和时间序列数据自动转化到MODFLOW程序的每个网格中，生成所需的数值模型；网格法需要为每个网格进行各类参数的赋值。但无论哪种方法，最终均生成MODFLOW需要的数值模型数据。TOUGH2的输入文件直接为数值模型数据。GMS生成的MODFLOW数值模型文件和TOUGH2输入文件见表 1。本次开发的转化程序是针对GMS生成出来的数值模型文件，而GMS各种版本生成出来的基本数值模型文件(表 1中列出的文件类型)相同，因此转化程序适用于所有GMS的版本。
| GMS/MODFLOW | 转换关系 | TOUGH2 | ||
| 文件后缀名 | 包含主要数据 |  |
关键字/文件 | 说明 |
| .dis | 分层高程和时间步长 | ROCKS | 岩性参数 | |
| .ba6 | 单元计算标示与初始水头 | PARAM | 模拟时间和步长 | |
| .lpf.mfs | 渗透系数、储水率或岩性编号 | ELEME、CONNE、MESH | 网格单元和连接 | |
| .wel | 井 | INDOM/INCON | 初始条件 | |
| .chd | 一类水头边界 | GENER | 井等源汇 | |
| .ghb | 通用水头边界 | FOFT | 设置单元格的结果输出 | |
| 2D/3D网格数据 | TIMES | 设置时间点的结果输出 | ||
| OUT | 结果文件 | |||
| 注：此表仅列出了本文所实现的转化文件。 | ||||
GMS/MODFLOW的空间剖分采用的是矩形网格，扩展名为.dis的文件包含着x和y方向的空间步长以及每层的顶底板高程数据。TOUGH2的网格文件MESH分为两部分：单元和连接。根据.ba6文件中的单元计算标示数据确定是否生成TOUGH2的单元和连接数据，然后读取.dis文件的数据计算体积、坐标、距离和面积等，读取.mfs文件赋岩性。具体转化过程见图 1。
|
| R1—R5.第1到第5行；C1—C4.第1到第4列；L1—L5.第1到第5层。dx1, dx2, dy1, dy2, dz1, dz2.沿x、y和z方向的连接中其第1和第2个网格到连接面的距离；dx, dy, dz.某网格在x、y和z方向上的长度。A1, A2.地层无起伏时的网格1和2相连的连接面面积。θ.连接面的法线与重力方向的夹角。 图 1 GMS/MODFLOW网格转化为TOUGH2网格MESH文件计算原理和储存格式 Figure 1 Method and storage format for conversion of the GMS/MODFLOW grid to the MESH file contains grid information |
|
|
值得注意的是，对于有起伏的地层，TOUGH2的以上处理方法与GMS/MODFLOW有不同之处。GMS/MODFLOW首先根据单元的厚度和渗透系数计算导水系数，然后采用几何平均法[18]计算相邻单元的界面导水系数；而TOUGH2在计算相邻单元的连接面积和渗透率时均采用算术平均法，等价于采用算术平均法计算导水系数。这个差异可能导致两者在计算水头时有一定的差异。
1.2.2 岩性参数转化GMS/MODFLOW和TOUGH2采用了不同的物理量表达岩性参数属性，如GMS/MODFLOW采用渗透系数K来刻画流体在地层中的流动能力，而在TOUGH2中采用的是渗透率k；GMS/MODFLOW中表达地层在水头变化时释水能力的储水率Ss需要转化为TOUGH2中的压缩系数α。转化关系为：
(1) 
(2) 式中：ρ为密度，此处近似取1 000 kg/m3；μ为动力黏滞系数，取1×10-3 kg/(m·s)；g为重力加速度，取9.8 m/s2；ξ为水的压缩系数，取4.5×10-10 Pa-1；由于孔隙度、岩体和水的压缩系数共同影响单相流体的流动，因此，孔隙度φ假定为0.2。
1.2.3 初始和边界条件转化GMS/MODFLOW中采用水头进行计算，而TOUGH2中采用的是压力，因此，水头需要转化为压力，其转化关系为
(3) 其中：patm为大气压力，取1.013 25×105 Pa；z为埋深。
初始条件转化时，由GMS/MODFLOW中的.ba6文件转化为TOUGH2中的INCON文件，需要考虑单元是否是有效的计算单元和边界单元。对于TOUGH2/EOS1中需要的温度，假定为25 ℃，气相饱和度设定为0。
边界条件有两种：一类水头边界和二类流量边界。TOUGH2中一类边界的刻画方法是给单元无穷大体积，因此，根据存放一类边界数据的.chd文件和通用水头边界的.ghb文件确定边界单元和其压力。二类边界GMS采用转化为井的方法刻画，下面介绍其转化。
1.2.4 源汇文件转化GMS/MODFLOW中的源汇有很多类型，例如井、蒸发、降雨补给、沟渠排泄等。实际上，不论何种形式的源汇，均可转化为井的注入或是抽出。出于此考虑，本文仅进行井的转化。GMS存放井的格式和转化后的TOUGH2文件格式如图 2所示。注(抽)序列时间点通过读入.dis文件应力期数据确定。在井的数据读入过程中应该注意，当读入井的数量为负数时表示与上一应力期的源汇相同。模拟时间和步长以及输出结果控制可根据应力期数据确定。完成转化后，根据模拟的需要修改或补充GMS不能创建的数据，然后进行计算，计算的结果可以转化为GMS的2D/3D格式数据进行显示。
|
| 图 2 GMS/MODFLOW井文件转化为TOUGH2源汇文件 Figure 2 Conversion of the well file in GMS/MODFLOW to the sink/source file in TOUGH2 |
|
|
转化的目的是借助GMS强大的前后处理能力建立复杂的TOUGH2模型。当然，我们可以进一步通过研究相同的模型来对比2个模拟器的模拟结果，分析模拟器的异同，从而增强模拟器的可信度。
2.1 水平多层变流量抽(注)水模拟对比模型为9层含(弱透)水层，边界为定水头边界，在中间位置有一抽(注)水井，滤网在第7和第9含水层位置(图 3)，抽(注)水量和模型相关参数根据某地区实际资料获得。
|
| 图 3 水平多层变流量抽(注)水模拟概念模型 Figure 3 Conceptual model for the water pumping (or injection) by variable mass flow in the horizontal multi-layer formation |
|
|
GMS/MODFLOW和由GMS/MODFLOW转化后TOUGH2的部分计算结果对比见图 4，可见2个软件的结果差别不大。
|
| 图 4 水平多层变流量抽(注)水模拟结果对比 Figure 4 Comparison of simulation results for the water pumping (or injection) by variable mass flow in multi-layer formation |
|
|
起伏非均质的模型根据GMS说明书第4个MODFLOW模型计算实例[19]修改而来。为了简化，修改初始水位为750 m，模拟设置为非稳定流，储水率为0.000 2 d-1。起伏地层GMS网格剖分见图 5，TOUGH2计算结果可以通过GMS导入散点并插值到2D网格进行显示，在此仅取模型第一层结果对比(图 6)；可见GMS/MODFLOW和TOUGH2计算结果总体趋势一致，存在一定的差异，可能是由于对起伏地层2个软件的处理方法不同所致。
|
| 图 5 基于GMS的起伏地层网格剖分 Figure 5 Mesh subdivision for the formation with ragged surface based on GMS |
|
|
|
| 图中显示为模型第一层在第151天时的水位。 图 6 起伏地层情况下GMS/MODFLOW和TOUGH2模拟结果对比 Figure 6 Comparison of simulation results from GMS/MODFLOW and TOUGH2 in the case of rolling formation |
|
|
1) GMS具有强大的前后处理功能，能够被用于建立复杂的三维模型，为TOUGH2快速建立数值模型。
2) GMS和TOUGH2虽然在模型的处理上有较大的不同，但是其模拟结果差异不大，说明两者均有很高的可靠性。
3) 虽然可以通过GMS为TOUGH2建立复杂的模型，但是现有的GMS数值模型的转化只针对于GMS/MODFLOW中常用的井、一类和通用水头边界模块，对于河流、沟渠、降雨补给和蒸发等模块的转化还待完善。
| [1] |
祝晓彬. 地下水模拟系统(GMS)软件[J].
水文地质工程地质, 2003, 30(5): 53-55.
Zhu Xiaobin. Groundwater Modeling System(GMS) Software[J]. Hydrogeology & Engineering Geology, 2003, 30(5): 53-55. |
| [2] | Owen S J, Jones N L, Holland J P. A Comprehensive Modeling Environment for the Simulation of Groundwater Flow and Transport[J]. Engineering with Computers, 1996, 12(3): 235-242. |
| [3] | Pruess K, Oldenburg C M, Moridis G J. TOUGH2 User's Guide[M]. 2nd ed. Berkeley: Lawrence Berkeley National Laboratory, 1999. |
| [4] |
郭亮亮, 张延军, 许天福, 等. 大庆徐家围子不同储层改造的干热岩潜力评估[J].
吉林大学学报(地球科学版), 2016, 46(2): 525-535.
Guo Liangliang, Zhang Yanjun, Xu Tianfu, et al. Evaluation of Hot Dry Rock Resource Potential Under Different Reservoir Conditions in Xujiaweizi Area, Daqing[J]. Journal of Jilin University (Earth Scicncc Edition), 2016, 46(2): 525-535. |
| [5] |
雷宏武, 金光荣, 李佳琦, 等. 松辽盆地增强型地热系统(EGS)地热能开发热-水动力耦合过程[J].
吉林大学学报(地球科学版), 2014, 44(5): 1633-1646.
Lei Hongwu, Jin Guangrong, Li Jiaqi, et al. Coupled Thermal-Hydrodynamic Processes for Geothermal Energy Exploitation in Enhanced Geothermla System at Songliao Basin, China[J]. Journal of Jilin University (Earth Science Edition), 2014, 44(5): 1633-1646. |
| [6] |
韩忠, 邵景力, 崔亚莉, 等. 基于MODFLOW的地下水流模型前处理优化[J].
吉林大学学报(地球科学版), 2014, 44(4): 1290-1296.
Han Zhong, Shao Jingli, Cui Yali, et al. Preprocessing Optimization of Groundwater Flow Model Based on MODFLOW[J]. Journal of Jilin University (Earth Science Edition), 2014, 44(4): 1290-1296. |
| [7] | Thunderhead Engineering. Petrasim User Manual [M]. 4th ed. New York:[s.l.], 2007. |
| [8] | Li Y, Niewiadomski M, Trujillo E, et al. Tougher: A User-Friendly Graphical Interface for TOUGHREACT[J]. Computers & Geosciences, 2011, 37(6): 775-782. |
| [9] | Florian Wellmann J, Croucher A, Regenauer-Lieb K. Python Scripting Libraries for Subsurface Fluid and Heat Flow Simulations with TOUGH2 and SHEMAT[J]. Computers & Geosciences, 2012, 43: 197-206. |
| [10] | Audigane P, Chiaberge C, Mathurin F, et al. A Workflow for Handling Heterogeneous 3D Models with the TOUGH2 Family of Codes: Applications to Numerical Modeling of CO2 Geological Storage[J]. Computers & Geosciences, 2011, 37(4): 610-620. |
| [11] | Borgia A, Cattaneo L, Marconi D, et al. Using a MODFLOW Grid, Generated with GMS, to Solve a Transport Problem with TOUGH2 in Complex Geological Environments: The Intertidal Deposits of the Venetian Lagoon[J]. Computers & Geosciences, 2011, 37(6): 783-790. |
| [12] |
杨艳林, 许天福, 李佳琦, 等. 应用TOUGH模拟二氧化碳地质储存过程的复杂地质体建模技术与实现[J].
吉林大学学报(地球科学版), 2014, 44(4): 525-535.
Yang Yanlin, Xu Tianfu, Li Jiaqi, et al. Complex Geological Body Modeling and Implementation of CO2 Geological Storage Simulation Using TOUGH[J]. Journal of Jilin University (Earth Science Edition), 2014, 44(4): 525-535. |
| [13] | Berry P, Bonduá S, Bortolotti V, et al. A GIS-Based Open Source Pre-Processor for Georesources Numerical Modeling[J]. Environmental Modelling & Software, 2014, 62: 52-64. |
| [14] | Bonduá S, Berry P, Bortolotti V, et al. TOUGH2 Viewer: A Post-Processing Tool for Interactive 3D Visualization of Locally Refined Unstructured Grids for TOUGH2[J]. Computers & Geosciences, 2012, 46: 107-118. |
| [15] | Avis J, Calder N. Walsh R. mView: A Powerful Pre-and Post-Processor for TOUGH2 [C]//Proceedings of TOUGH Symposium. Berkeley: Lawrence Berkeley National Laboratory, 2012: 1-11. |
| [16] | Hu L, Zhang K, Cao X, et al. IGMESH: A Convenient Irregular-Grid-Based Pre-and Post-Processing Tool for TOUGH2 Simulator[J]. Computers & Geosciences, 2016, 95: 11-17. |
| [17] | Harbaugh A W, MODFLOW-2005: The U S Geological Survey Modular Ground-Water Model:The Ground-Water Flow Process, in Techniques and Methods [R]. Reston: U S Geological Survey, 2005. |
| [18] | McDonald M G, Harbaugh A W. A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model [R]//Techniques of Water-Resources Investigations of U S Geological Survey. Reston: U S Geological Survey, 2005. |
| [19] | Environmental Modeling Research Laboratory. Groundwater Modeling System Tutorials: Vol 2: Modflow-Conceptual Modeling Approach[R]. Provo: Aquaveo, 2004. |