松辽盆地增强型地热系统（EGS）地热能开发热－水动力耦合过程

雷宏武，金光荣，李佳琦，石 岩，冯 波

吉林大学地下水资源与环境教育部重点实验室，长春 130021

0 引言

干热岩（hot dry rock，HDR）作为地热资源的一种，由于分布广泛、清洁和热储量巨大，被认为是21世纪最有潜力的资源［1］。针对内部不存在流体或仅有少量地下流体且渗透率极低的高温岩体（干热岩），增强型地热系统（enhanced geothermal system，EGS）或工程型地热系统（engineered geothermal system）作为一种新型的地热发电技术，并不需要天然热液资源的存在，可通过水力压裂的方式从干热岩中提取地热资源［2］。在地热资源的提取过程中，温度场、压力场、力学场和化学场互相影响，并随着时空发生复杂的变化。采用传统的方法无法全面了解这种复杂耦合过程的互相影响，而数值模拟作为地质资源和环境研究的新方法［3］，在耦合过程的研究中发挥着越来越大的作用。Yano和Ishido［4］采用数值模拟方法研究了深部地热在超临界条件下的生产行为，表明压力在超临界区域变化复杂，其主要受流体动力黏滞性和压缩能力的非线性特征控制。FRACTure和CHEMTOUGH被耦合形成FRACHEM，用于研究了Soultz深部地热的多场耦合问题［5－7］。Xu等［8］采 用 TOUGHREACT研究了Soultz深部地热系统注入不同成分水的地球化学演化特征，以及其对注入能力的影响。Pruess［9－10］采用 TOUGH2对比分析了 CO2和水作为载热介质时的地热系统中的压力和温度变化特征，以及热提取效率差别。

笔者以松辽盆地徐家围子为研究对象，采用数值模拟方法研究了其深部地热开发过程中温度和压力的变化特征，并分析了模型存在的各种不确定性，以期为下一步工程实施提供理论依据。

1 松辽盆地地质概况

松辽盆地位于中国东北部，总面积约26万km2，是一个大中型中新生代陆相沉积的、呈北北东向展布的菱形盆地。盆地主要有6个一级构造单元［11］（图1）：中央凹陷区、东北隆起区、东南隆起区、西南隆起区、西部斜坡区、北部倾没区。盆地主要沉积了中、新生代的沉积岩系，最大沉积厚度大于10 km，地层从下而上依次为上侏罗统、白垩系、第三系和第四系。盆地深层是指泉二段以下地层，自上之下由泉头组、登娄库组、营城组、沙河子组和火石岭组组成［12］。

松辽盆地现今平均地温梯度是3.7℃/100m，是中国大型盆地中地温梯度比较高的地区之一［13］。徐家围子断陷位于中央凹陷区东北部，主要由徐西断裂、宋西断裂和徐东断裂分割成2个较大的断裂区（图2），是松辽盆地深部重要的油气储层，并且其靠近相对较高的地温梯度带，是优良的深部地热开采靶区之一。

宋西断裂和徐东断裂构成的断裂区沉积厚度较大，约4.5km，根据 MIT报告［1］对增强地热系统的定义，适用于深部地热开采的深度为3.0～10.0 km，因此，本区营城组中下部、沙河子组和火石岭组将是深部地热开采的备选层位。营城组为火山岩－陆源碎屑岩沉积建造，下部主要由安山玄武岩、火山角砾岩、砂岩及砂砾岩组成，上部为酸性火山岩、火山碎屑岩及砂岩和粉砂岩；沙河子组以泥岩为主，夹有砂岩、粉砂岩及少量凝灰岩；火石岭组主要以凝灰质角砾岩、凝灰岩、安山岩、玄武岩及凝灰质砾岩为主。

图2 徐家围子断陷地质剖面概化图（据文献［12］修改）Fig.2 Geologic profile across the Xujiaweizi depression in E－W direction（modified after reference［12］）

2 模型建立

2.1 模拟工具

用于地热开采的数值模拟软件可分为2类：干热岩和水热型地热系统模拟器。干热岩中典型的模拟 软 件 有 FRACTure［14］、 GEOTH3D［15］、FRACSIM－3D［16］和 Geocrack2D［17］等，水热型中常用的有 TOUGH2［18］、STOMP［19］、和 FEHM［20］等。TOUGH系列由于其能够刻画多个过程和可信的模拟结果［21］，被广泛用于与地热有关的研究中［8，22－23］，其灵活的网格剖分方法和多重连续介质模型［24］（multiple interacting continua，MINC）使得其不仅适用于水热型地热系统，而且也能够用于采用水作为载热流体的孔隙和裂隙共同存在的干热岩地热系统中。因此，本次研究拟采用TOUGH系列软件进行。

TOUGH2为TOUGH系列的主要成员，其基于组分质量和能量守恒求解多相流动过程中压力、温度、饱和度和组分质量分数的时空分布，采用的空间离散方法是积分有限差，时间离散是隐式差分，非线性方程求解是Newton－Raphson迭代［18］。TOUGH2采用模块化设计，EOS1模块能够刻画纯水、水蒸气及两相共存系统，温度最高可达647.3 K，其两相饱和线、密度、黏度、热焓等根据国际公式化委员会提供的模型计算。因此，TOUGH2/EOS1可用于松辽盆地深部地热开采的模拟工作，其主要特征见文献［18］。生产井和注入井的刻画采用生产指数模型［25］。

2.2 模型概化

深部热储层孔隙度和渗透率一般都较低，往往需要进行人工压裂形成新的裂隙来增强系统的流动能力。裂隙系统提供了流体快速流动的通道，保证了热量可以循环不断的提取。徐家围子断陷区深层为沉积岩和火山岩混杂，具有一定的孔隙度和渗透率，可以认为是连续的孔隙系统。因此，通过人工改造以后的徐家围子深部地热系统是一个孔隙－裂隙复合的复杂系统。

本次研究将徐家围子深部地热系统在垂向上概化为孔隙和裂隙2个互相关联的统一系统。假设压裂范围垂向上为50m，采用MINC刻画。热提取过程中由于上下孔隙系统的存在，注入井附近流体会向孔隙系统漏失，生产井附近流体会向裂隙系统流入，同时孔隙系统的热会向裂隙系统进行传导。为了刻画这种影响，裂隙系统上下各取500m的孔隙介质。平面上为了简化边界的影响，采用“五点法”布井方式［9－10］（图3a）。由于对称性，模型只需取1/8进行计算。平面上模型网格剖分间距为70.7m×70.7m，共剖分36个网格；垂向上网格间距为50.0 m，共21层。裂隙系统网格为180个，孔隙系统为720个，模型共900个网格（图3b）。

徐家围子断陷深部可供地热开采的层位有3层，为了分析开采过程中不同层位热和水动力的异同，笔者建立了3个典型位置的模型（埋深为3 250、3 750、4 100m，分别以 Case1、Case2和 Case3代表），其对应的中间压裂层的温度分别为131、147和160℃，初始压力符合静水压力，整个模型的初始温度按照T＝111＋0.035（D－2 500）（D＞2 500m）给定。注入井和生产井均采用定压力方式，分别高于和低于初始压力2MPa，模拟开采时间为36a。模型具体参数见表1。

图3 “五点法”地热开采概念模型和网格剖分Fig.3 Schematic representation of the five－spot well pattern geothermal development with computational grid

表1 “五点法”开采热储层模型基本物理参数Table 1 Basic property of geothermal storage formation in five－spot well pattern

3 模拟结果及分析

3.1 压力场时空变化特征

从压力时空化特征（图4）可以看到：裂隙系统的压力相对初始压力大部分是降低的，生产井附近开采引起的压力降低较大，而注入井由于持续注入压力降低较小；整个裂隙系统的压力分布明显受生产井的控制，但随着时间的增加，注入的影响逐渐增大，整个系统压力逐渐上升。这种变化特征主要原因是水的移动性（ρ/μ，其中：ρ和μ分别为流体密度和黏度）在注入和生产区域的不同，注入区域由于冷水的注入使得水的温度降低，导致其移动性比生产区域的降低了约80%，从而使注入区域的压力梯度明显大于生产区域。随着热能逐渐被开采，生产井附近温度降低，水的移动性的差别减小，注入的影响增大。

从图4不同埋深位置的结果可以看到，埋深越大，整个系统的压力降低程度越大，其主要原因是埋深大，温度高，水的移动性增强，定压产出条件下的水的流出质量大，从而使得压力降低大。

3.2 温度场时空变化特征

注入的冷水在由注入井到生产井流动的过程中，与其接触的介质进行对流传导换热，由于裂隙介质体积比例远远小于孔隙介质，因此，换热过程主要发生在流体和孔隙介质之间。孔隙介质的低孔低渗特征又使换热方式以传导为主。

温度变化由注入井向生产井是逐渐减小的；随着时间的增加，由生产井到注入井的地热逐渐提取完毕（图5）。由不同埋深的模拟结果（图5）对比可以看到：随埋深增大，温度变化增大，主要原因是埋深大的热储层初始温度较高；另外其温度变化影响范围要大于埋深小的，主要由于温度高加快了水的流动。

从图6中孔隙基质温度的变化特征可以看到：与流体温度变化特征相似，由生产井到注入井温度逐渐降低；并且孔隙基质中温度分布存在梯度，随着与裂隙距离的增加，温度降幅减小，反映了孔隙基质中由内向外的传热过程。孔隙基质中的这种温度梯度在整个温度变化的前缘较大，而在后缘逐渐变小，反映了后缘部分孔隙基质中的热逐渐被提取完毕。埋深大的温度变化大，同样是由于孔隙基质本身的初始温度较高，加快了水的流动速度。裂隙系统上下未压裂处理的孔隙基质的温度变化不大，主要变化集中在压裂区域上下附近，并且上下变化量很相近，说明未压裂和压裂区域之间的热交换以传导为主，对流传热可以忽略。

图4 不同埋深位置地热开采过程压力变化时空分布Fig.4 Spatial and temporal distribution of pressure at different depth levels

图5 不同埋深位置地热开采裂隙介质中温度变化时空分布Fig.5 Spatial and temporal distribution of temperature in fractured media at different depth levels

图7显示了生产井和注入井连线（最短流动路径）的温度和压力变化。可以看到：注入井区域压力上升范围有限，大约50m，这主要是受注入井低温的影响；注入井附近压力梯度明显大于生产井区域；不同埋深情况的压力变化并不是很大，但温度变化较大。

从图8可以看出：生产井附近压力由于生产井压力降低而迅速降低并稳定，在开采后期整个系统的温度降低，温度引起的流动阻力增大，从而压力略有上升；埋深大，初始温度和压力均较高，从而在开发过程中其也相对较高。

3.3 热提取速率

从图9可以看到：对于埋深3 250m（Case1），由于开采过程中温度逐渐降低，注入质量速率和生产质量速率由开始的70kg/s左右逐渐减小到53kg/s左右，并且注入质量速率稍大于生产质量速率；埋深越大，注入质量速率和生产质量速率均增大，相对埋深3 250m而言，埋深3 750m（Case2）和4 100m（Case3）的速率分别增大了2%～4%和4%～6%。埋深3 250m的热提取速率由开始的40MW左右减小到15MW左右（图10）。与质量速率变化一样，热提取速率也随着埋深增大而增加，埋深3 750 m和4 100m的热提取速率分别增加了11%～17%和19%～29%，明显大于质量速率的增幅，主要原因是埋深大的流体温度高，本身携带的能量相对较大。

因此，从热提取速率的角度来说，松辽盆地徐家围子3个埋深水平中埋深4 100m水平明显要优于其他埋深水平，同时也说明了埋深越大，温度越高，热提取速率也高。但深度增加，也必然增加了钻井成本。

4 参数不确定分析

地质条件的复杂性使得参数存在很多的不确定性，从而导致模拟结果的不确定性。为了评价这种不确定性，笔者进行了孔隙基质渗透率和孔隙度、裂隙介质渗透率和裂隙度、岩石导热系数、井径、注入/生产压力、注入温度等因素的不确定性分析，具体模拟设置见表2。

图6 不同埋深位置地热开采注入井－生产井之间孔隙介质中温度变化时空分布Fig.6 Spatial and temporal distribution of temperature in porous media at different depth levels

图7 注入井－生产井之间压力和温度变化量分布（36.5a）Fig.7 Change of pressure and temperature profiles along a line from injection to production well after a simulation of 36.5years

图8 生产井所在网格压力和温度时间变化特征Fig.8 Evolution of pressure and temperature of production well

图9 注入井和生产井质量速率Fig.9 Mass rate of injection and production well

图10 热提取效率时间变化特征Fig.10 Rate of heat extraction and cumulative heat produced for five－spot problem

4.1 孔隙介质渗透率和孔隙度

从图11和图12可以看出：孔隙基质的渗透率和孔隙度对压力和温度的影响很小，整体上在靠近生产井区域压力随着孔隙基质渗透率和孔隙度的增大稍微有所降低，温度随孔隙基质渗透率和孔隙度的增大稍微有所升高；质量和热提取速率整体上是随着孔隙基质渗透率和孔隙度的增大而升高的。

4.2 裂隙介质渗透率和裂隙度

从图13可以看出：裂隙介质渗透率对温度、压力、质量和热提取速率的影响很大；裂隙介质渗透率逐渐增大时，注入井逐渐占主导地位，压力增大，温度降低；质量和热提取速率随着裂隙介质渗透率的增大而增大，当裂隙渗透率增加到10－12数量级时，热提取速率变化剧烈，由开始的160MW到开采15a时的10MW，热储存中的热几乎被开采完毕。因此，裂隙介质渗透率直接决定了压力和温度的时空分布特征，以及质量和热提取速率。

表2 不同条件的数值模拟设置Table 2 Numerical simulation runs for different conditions

从图14可以看出：裂隙介质裂隙度对温度、压力、质量和热提取速率没有影响。这是因为裂隙介质所占比例很小，裂隙度的变化带来的影响可以忽略，但是裂隙度直接影响流体的实际流速，从而可以改变整个系统的循环周期。

4.3 热传导系数

从图15可以看出：压力几乎不受岩石的热传导系数的影响；温度降幅随着热传导系数增加而减小；而质量和热提取效率相应地增加。这是因为热传导系数增加使得基质中的热更快地传导到裂隙介质中来，冷水注入引起的温度降低减小，加快了流体的流动过程，从而提高了质量和热提取速率。但因岩石的热传导系数间差别较小，这种影响并不明显。

图11 孔隙介质渗透率对温度、压力和质量及热提取效率的影响Fig.11 Effects of permeability of porous media on temperature，pressure，and mass and heat extraction rate

4.4 井径

从图16可以看出：随着井径的增加，井与地层的接触面积增大，使得注入速率增大，导致压力降低程度减小，温度在靠近生产井区域降低程度增大。在开采过程中，随着井径的增大，质量速率增大，而热提取速率在开采早期增大，开采后期这种差别逐渐减小。这主要是因为井径大的开采后期时系统温度降低快，故单位质量的流体携带的热量变小。

4.5 注入井/生产井压力和温度

图12 孔隙基质孔隙度对温度、压力和热提取效率的影响Fig.12 Effects of porosity of porous media on temperature，pressure，and mass and heat extraction rate

从图17可以看出：注入井/生产井压力大，在注入井区域（x＝0.0m）压力增大得多，在生产井区域（x＝700.0m）压力降低得多，生产井区域的变化幅度明显大于注入区域，生产井控制整个系统的压力，同时开采质量和热提取速率也高，温度降低程度大；注入/生产压力大的，随着热被迅速提取，温度降低，其热提取速率迅速变小。

从图18a可以看出：注入温度的增加可以增强注入区域水的流动性，提高注入井的影响范围，从而使得整个系统的压力降低较小，大大提高质量流动速率；但是注入的能量明显增加，最终的能量净提取效率并不一定增加。图18b显示随着注入温度的增加，质量速率增加，但热提取速率整体减小。温度分布在注入井附近注入温度越大，降低越小，而在生产井附近由于注入温度高的质量速率大，从而其温度降低得快。

图13 裂隙介质渗透率对温度、压力和热提取效率的影响Fig.13 Effects of permeability of fractured media on temperature，pressure，and mass and heat extraction rate

4.6 基质

从图19可以看出：不考虑压裂区域外孔隙介质的影响，温度降低明显增大，导致了质量流速和热提取速率的降低，压力变化幅度减小。可见周围孔隙介质对压裂区域传导的热量还是占有一定比例的。因此，在深部地热开采模拟过程中应该考虑上下基质传热过程。

图14 裂隙介质裂隙度对温度、压力和热提取效率的影响Fig.14 Effects of porosity of fractured media on temperature，pressure，and mass and heat extraction rate

通过参数的不确定分析可知，裂隙系统的渗透率、注入井/生产井压力和注入温度、井径对深部地热开采过程中的压力和温度影响较大，从而影响到热的提取效率，应重视这些因素对松辽盆地地热开采过程的影响。而孔隙介质的渗透率和孔隙度、裂隙介质的孔隙度和岩石的热传导系数的影响并不明显。

5 结论与展望

1）松辽盆地徐家围子断陷区地温梯度较高，埋深4 100m地层温度达到160℃，是良好的热储层；其岩性主要为砂岩，具有一定的孔隙度和渗透率（孔隙度在0.1左右，渗透率在10－17m2数量级上），但仍然需要采用水力压裂来人工增加渗透性，以达到规模化开发的目的。

图15 岩石热传导系数对温度、压力和热提取效率的影响Fig.15 Effects of rock thermal conductivity on temperature，pressure，and mass and heat extraction rate

2）在采用定压力开采地热的情况下，注入井和生产井区域温度的差异使得水的流动性差别较大，生产井控制整个区域的压力分布，压力梯度在注入井区域较大，随着开发的进行，注入井的影响逐渐增大；裂隙介质中的流体温度由注入井到生产井逐渐升高，并随着开发的进行而降低，逐渐向生产井发展，其中的孔隙基质温度变化特征与流体相似，并且在温度变化前缘存在明显的梯度，即离裂隙越近，温度降低越大。

3）不同埋深位置的模拟显示，埋深大的温度相对较高，水的流动性较强，质量和热提取速率较高，压力和温度变化均较大。

图16 井径对温度、压力和热提取效率的影响Fig.16 Effects of well radius on temperature，pressure，and mass and heat extraction rate

4）裂隙系统的渗透率、注入井/生产井压力和注入温度、井径对深部地热开采过程中的压力和温度影响较大，从而影响到热的提取速率，而孔隙基质的渗透率和孔隙度、裂隙介质的孔隙度和岩石的热传导系数的影响并不明显。

本次研究结论是在特定参数条件下获得的，并且只考虑了压力场和温度场的耦合效应，没有综合考虑应力场和化学场的影响，而实际工程条件下这四个场的互相影响是客观存在、复杂和重要的，因此，还需在以后的研究中进一步深入。

图17 注入井/生产井压力对温度、压力和热提取效率的影响Fig.17 Effects of injection/production well pressure on temperature，pressure， and mass and heat extraction rate

图18 注入温度对温度、压力和热提取效率的影响Fig.18 Effects of injection temperature on temperature，pressure，and mass and heat extraction rate

图19 压裂区域上下基质对温度、压力和热提取效率的影响Fig.19 Effects of surrounding matrix on temperature，pressure，and mass and heat extraction rate

（References）：

［1］Massachusetts Institute of Technology.The Future of Geothermal Energy：Impact of Enhanced Geothermal Systems（EGS）on the United States in the 21st Century［R］.Cambridge：MIT，2006.

［2］Lund J W.Characteristics，Development and Utilization of Geothermal Resources［J］.Geo－Heat Center Bulletin，2007，28（2）：1－9.

［3］许天福，金光荣，岳高凡，等.地下多组分反应溶质运移数值模拟：地质资源和环境研究的新方法［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.

［4］Yano Y，Ishido T.Numerical Investigation of Production Behavior of Deep Geothermal Reservoirs at Super－Critical Conditions［J］.Geothermics，1998，27（5/6）：705－721.

［5］Bachler D，Kohl T.Coupled Thermal－Hydraulic－Chemical Modelling of Enhanced Geothermal Systems［J］.Geophysical Journal International，2005，161（2）：533－548.

［6］Rabemanana V，Durst P，Bachler D，et al.Geochemical Modelling of the Soultz－Sous－Forets Hot Fractured Rock System：Comparison of Two Reservoirs at 3.8 and 5km Depth［J］.Geothemics，2003，32（4/5/6）：645－653.

［7］Andre L，Rabemanana V，Vuataz F D.Influence of Water－Rock Interactions on Fracture Permeability of the Deep Reservoir at Soultz－Sous－Forets，France［J］.Geothemics，2006，35（5/6）：507－531.

［8］Xu T F，Sonnenthal E，Spycher N，et al.TOUGHREACT－A Simulation Program for Non－Isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media：Applications to Geothermal Injectivity and CO2Geological Sequestration［J］.Computers ＆ Geosciences，2006，32（2）：145－165.

［9］Pruess K.Enhanced Geothermal Systems（EGS）U－sing CO2as Working Fluid：A Novel Approach for Generating Renewable Energy with Simultaneous Sequestration of Carbon［J］.Geothermics，2006，35（4）：351－367.

［10］Pruess K.On Production Behavior of Enhanced Geo－thermal Systems with CO2as Working Fluid［J］.Energy Conversion and Management，2008，49（6）：1446－1454.

［11］高瑞琪，蔡希源.松辽盆地油气田形成条件与分布规律［M］.北京：石油工业出版社，1997.Gao Ruiqi，Cai Xiyuan.Field Formation Conditions and Distribution Rules in Songliao Basin［M］.Beijing：Petroleum Industry Press，1997.

［12］周庆华，冯子辉，门广田.松辽盆地北部徐家围子断陷现今地温特征及其与天然气生成关系研究［J］.中国科学：D辑：地球科学，2007，37（增刊II）：177－188.Zhou Qinghua，Feng Yuhui，Men Guangtian.The Present Temperature Characteristics and Relationship Between Natural Gas Generation of Xujiaweizi Fault in Northern Songliao Basin［J］.Science in China：Series D：Earth Science，37（Sup.II）：177－188.

［13］任战利.中国北方沉积盆地构造热演化史研究［M］.北京：石油工业出版社，1999.Ren Zhanli.Study on the Tectonic－Thermal History of Sediment Basin in Northern China［M］.Beijing：Petroleun Industry Press，1999.

［14］Kohl T，Hopkirk R J.“FRACTure”：A Simulation Code for Forced Fluid Flow and Transport in Fracture，Porous Rock［J］.Geothemics，1995，24（3）：333－343.

［15］Yamamoto T，Kitano K，Fujimitsu Y，et al.Application of Simulation Code，GEOTH3D，on the Ogachi HDR Site［C］//Procedings of the 22nd Annual Workshop on Geothermal Reservoir Engineering.Stanford：Stanford University，1997.

［16］Jing Z，Willis－Richards J，Watanabe K，et al.A Three－Dimensional Stochastic Rock mechanics Model of Engineered Geothermal Systems in Fractured Crystalline Rock［J］.Journal of Geophysical Research，2000，105（B10）：23669－23679.

［17］Swenson D.User’s Manual for Geocrack2D：A Coupled Fluid Flow/Heat Transfer/Rock Deformation Program for Analysis of Fluid Flow in Jointed Rock［R］.Manhattan：Kansas State University，1997.

［18］Pruess K，Curt O，George M.TOUGH2USER’S GUIDE，VERSION 2.0［R］.Berkeley：University of California，1999.

［19］White M D，Oostrom M.STOMP－Subsurface Transport over Multiple Phases，Version 4.0，User’s Guide［R］.Richland：Pacific Northwest National Laboratory，2006.

［20］Zyvoloski G A，Robinson B A，Dash Z V，et al.User’s Manual for the FEHM Application：A Finite－Element Heat－and Mass－Transfer Code［R］.Los Alamos：Los Alamos National Laboratory，1997.

［21］Pruess K，García J，Kovscek T，et al.Code Intercomparison Builds Confidence in Numerical Simulation Models for Geological Disposal of CO2［J］.Energy，2004，29：1431－1444.

［22］Croucher A E，O’Sullivan M J.Application of the Computer Code TOUGH2to the Simulation of Supercritical Conditions in Geothermal Systems［J］.Geothermics，2008，37（6）：622－634.

［23］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（4）：429－442.

［24］Pruess K.GMINC－A Mesh Generator for Flow Simulations in Fractured Reservoirs［R］.Berkeley：University of California，1983.

［25］Coats K H.Geothermal Reservoir Modeling［C］//The 52nd Annual Fall Technical Conference and Exhibition of the SPE.Denver：Society of Petroleum Engineers，1977.