增强型地热系统研究综述

郭剑，陈继良，曹文炅，蒋方明

(1.中国科学院广州能源研究所，广州市510640;2.中国科学院可再生能源重点实验室，广州市510640)

0 引言

目前，我国的电力结构是以火力发电为主，在提供能源供给的同时也造成了较大的环境问题，而且煤炭、石油等化石能源正在日益耗竭，因此寻找环境友好的新型能源是我国能源领域的重要任务。2013年1月，国家能源局等部门联合发布了《关于促进地热能开发利用的指导意见》，指出我国地热能资源储量巨大，明确提出了未来地热能发展目标:到2015年，地热发电装机容量达到100 MW，地热能利用量达到2 000万 t标准煤;到2020年，地热能利用当量达5 000 万 t标准煤［1］。

地热能具有稳定、连续、利用系数高等优势，将会成为未来能源中的重要部分。但是贮存在地下3 km以内、温度＜150℃的中、低温地热资源仅占地球地热的10%，以水热型为主，用于发电则效率相对较低，而且热田规模通常较小，难于维系长时间稳定产电要求。旨在开发和利用地底深层3～10 km、以干热岩(hot dry rock，HDR)热能为主的增强型地热系统(enhanced geothermal systems，EGS)正在逐渐引起人们的广泛关注［2］。

鉴于EGS巨大的开发应用前景，世界各国正在展开广泛的研究，一些国家先后建立了野外试验场或示范电站，例如美国的 Fenton EGS、法国的 Soultz EGS、英国的 Rosemanowes EGS、日本的 Hijiori EGS、澳大利亚的Cooper EGS等。这些国家经过近40年的野外试验研究，已经在钻井探测、水力压裂、人造热储和采热循环等方面获得了技术成就。而我国的EGS研究还处于起步阶段，目前尚未开展野外试验研究，只进行了一些资源评价［3-4］、前景分析［5-11］方面的研究工作。本文通过文献综述对国外野外试验场的研究经验以及国内外数值模拟的研究方法进行整理和介绍，以期为我国的EGS研究和电站建设提供参考。

1 EGS原理

通过人为的手段如水力激发(hydraulic stimulation)，在低渗透性结晶质HDR内建设人工热储，然后灌注冷流体介质、置换HDR的热能，将热流体开采出来后用于地面发电，如图1所示。

图1 增强型地热系统建设过程示意Fig.1 Procedure of EGS construction

EGS产热温度通常可达150～350℃，用于发电具有高能量转换效率，而且干热岩的地热资源储量也巨大。麻省理工学院在2006年的一份报告［12］指出，美国干热岩地热资源总储量为14×106EJ(不包括黄石公园地区)，如果保守估计2%的有效开发利用率，基于2005年美国总能耗100 EJ，干热岩地热资源可满足美国2 800年的能耗需求。该报告同时预测，只要投入足够的经费开展基础研究和应用技术研发，到2050年，美国干热岩地热电站将实现产能100 000 MW以上。中国科学院地质与地球物理研究所汪集旸院士等的研究数据［3］显示，我国大陆地区干热岩地热储量为21×106EJ，而2010年我国总能耗为95.2 EJ，如果同样假定2%的开发利用率，干热岩地热资源可以满足我国4 400年左右的能耗所需。

2 野外试验场

EGS野外试验场建设包含钻井、水力激发、循环采热等多个阶段，通过井下探测、微震监测和水力循环参数等来研究相关问题。在众多的野外试验场中，以Soultz EGS的研究经验和取得成果最为丰富。本文以Soultz试验场的研究为例，依次从井探研究、水力激发、微震监测和循环测试这些方面对相关研究进行综述。

2.1 井探研究

井探数据主要包括切削样本(chip samples)、岩石碎片(rock fragments)、光谱射线(spectral gamma ray)和超声井像(ultrasonic borehole imagery，UBI)等。除了对切削样本和岩石碎片的直接观察获取的矿物组成、岩石纹理、水热填充物等初步信息，对其进行X射线衍射分析可以比较准确地区别蚀变裂隙区域和完整岩体区域。井下γ射线光谱分析用于确定岩相主要分布情况，包括种类变化、标准岩基区域和水热蚀变裂隙区域。用分层上升归类(hierarchy ascending classification，HAC)方法对γ射线光谱数据处理获得连续的井下剖面图，然后去除噪声获得光滑的HAC，与原始的HAC比较可以显示裂隙的存在迹象。裂隙的高精度位置、倾斜方向和倾角信息可从超声井像的波形和传输时间的分析获得，例如根据振幅和传输时间确定裂隙张开与否，根据正弦连续性确定裂隙的大小尺度。超声井像数据已成为井下裂隙探测的主要依据。

2.1.1 应力状态

根据井探数据对小尺度裂隙的研究，可以按照裂隙的初始状态将其划分为:(1)自然裂隙(natural fracture)，在钻井之前即已存在，与花岗岩岩体的地址构造时代有关;(2)诱发裂隙(induced fracture)，在钻井过程中形成，与当前应力场状态有关［13］。通过诱发裂隙特征推测应力场状态是应力研究的重要方式。Cornet和 Jones［14］经过研究，确定了地下 1.5～3.5 km的最大横向应力方向为N170°E±15°，该结论获得了大部分研究者的认可［15-16］。近来，Valley和Evans［17］根据各井最新的数据，将地下1.5～5 km 的最大横向应力方向确定为N169°E±14°，仅存在微小的差异，但是在应力大小与地理深度的线性关系上有较大不同。

从整个断层特征上看，Soultz地区存在2个主要的节理组，走向分别为 N10°E和 N170°E，倾向分别为65°W、70°E，节理的平均密度高达3.2个/m，反映该地区花岗岩中存在普遍断裂的现象［18-20］。

2.1.2 裂隙特征

图2 Soultz各井下裂隙方向岩相剖面图Fig.2 Major petrographical crystalline sections in wells of Soultz

根据Genter等［21］、Valley 等［22］、Dezayes 等［23-24］以及其他一些学者［25-27］的研究，各井地下裂隙结构特征基本已获认知。图2显示了各井下花岗岩层中裂隙方向信息［24］。可以看出，花岗岩中裂隙主要方向以较强的东向和西向的趋势从N160°E变化到NS。Valley等［22］用统计学方法将所有裂隙分为7组，其中60%的裂隙属于N-S走向的倾东和倾西这2组，此外还有NW-SE走向和NE-SW走向的这2个次垂直组，以上4组涵盖了花岗岩中已发现的95%的裂隙。对于GPK3和GPK4井，在实际垂直深度(true vertical depth，TVD)为1 420～2 700 m，裂隙组主要向东倾斜;在井中部TVD为2 700～4 800 m，主要是向东和向西倾斜的共轭裂隙组;在井底TVD为4 800～5 000 m，裂隙组主要是向西倾斜。

Dezayes等［23-24］对GPK3和GPK4井的裂隙密度随着深度变化进行了研究。GPK3井可以按照不同裂隙密度值分为9段，其中最大密度值为1.4 fract/m，出现在3 000～3 100 m以及4 800～4 900 m，仅占很小一部分，大部分深度以0.3、0.4 fract/m的小密度值为主。对于GPK4井，出现了密度值高达2.86、1.76 fract/m的深度段，也仅占井深的微小部分，而在3 500～4 750 m内，裂隙密度以最低值0.23 fract/m维持不变。

由于技术方法的局限，超声井像只能探测到裂隙的位置、倾角和倾向等信息，而无法探测裂隙的种类、水热填充物的性质以及裂隙的长度等信息，另外在裂隙的厚度探测方面也存在较大的困难。Genter等［21］对参考井ESP1进行了连续的岩心分析，对其地下3 000多个裂隙和一些裂隙区域进行了表征，尤其是对岩心裂隙的厚度进行了研究。其结论是裂隙的厚度按照幂函数分布，大小为0.1～250 mm。Genter等［21］把这些裂隙按照形成机理分为剪切断裂的延伸裂隙(Mode I)和纯粹的剪切断裂裂隙(Mode II)，厚度为0.1～10 mm的小裂隙主要属于前者，而其他大部分裂隙都属于后者。

2.1.3 裂隙区域

目前，在Soultz地区的6口井下一共确定了39个裂隙区域，如图3所示［21，24，28］。可以将其划分为3个级别，每个级别反映了不同的渗流能力。一级(Level I)涉及到主要裂隙区域，在钻井过程中已探测到的并伴随着严重的泥浆损失现象，在水力激发之前即存在一定的渗流能力;二级(Level II)裂隙区域在水力激发时显示了20%的水流损失，其特征为1条以上的厚裂隙伴随着一些水热蚀变的圈痕;三级(Level III)裂隙区域在水力激发时水流损失在20%以下，基本没有水热蚀变的圈痕。

在地下约1 200 m处有1个Level I级的大断层裂隙区域，分为3个子区域，分别和EPS1、GPK1、4550这3口井相交，GPK1和4550井的钻井过程泥浆损失主要就发生在此处。对于GPK1井，地下约1 800 m处还有1个Level I的主裂隙区域，在钻井过程中也发生过卤水自发渗流。对于GPK2井只在测量深度(measured depth，MD)2 120 m处发现1个Level I级裂隙区域，在钻井过程中所有的泥浆损失都发生在这里，并且在后来的水力激发过程中95%的水流从这里渗入。在MD为2 120 m以下，各井都没有Level I裂隙区域，只在GPK3井底部MD为4 770 m处例外，GPK3井70%的渗流对应此区域［20］。

图3 Soultz各井下裂隙区域综合图Fig.3 Synthesis of all fracture zones along Soultz wells

总体上看，这39个裂隙区域可以按照不同的深度范围划归为3个群集:TVD为1 800～2 000 m的Cluster I，TVD 为3 000～3 400 m 的Cluster II和TVD为4 500～5 000 m的Cluster III。这种划分除了对聚合特征不明显的GPK4井外，对其他各井都比较吻合。Genter等［29-30］研究倾向认为花岗岩中的这3个群集反映了1个较大的断层，与在沉积层中发现的断层可能等同。Cluster I位于斑岩花岗岩层中，包含有Level I的具有一定渗透率的主要裂隙区域。Cluster II未包含Level I的裂隙区域，位于有蚀变断裂的花岗岩层中。该岩层由于普遍和强烈的蚀变，具有大量中小尺度的裂隙，并且导致岩体弱化，从而易于受水力激发提高渗透率［31］。1997年的水力循环试验表明GPK1和GPK2井通过该区域形成了良好的连通［32］。Cluster III所在的岩层也具有一定低程度的普遍蚀变，该区域也是目前热储层底部主要的水力通路。

裂隙区域由众多裂隙复合构成。在某些情况下(例如，当裂隙无法完全看做一个平面时)，裂隙的方向就变得难以计算。因此通过裂隙方向的合成来确定裂隙区域的方向也就比较困难。通常在裂隙区域中有一些特殊的角砾化和微角砾化的裂隙轨迹，在一定程度上这些角砾化和微角砾化区域的接界可以假设为裂隙区域的方向。通过这种假设方法，一些关于裂隙区域方向的结论也被一些研究者给出［23-24］，但是他们也都强调了不准确性。

一般研究把裂隙区域的群集现象归因于莱茵地堑构造中发生的大尺度常规断裂［30］。但是Dezayes等［23］指出，裂隙区域的方向以及中尺度裂隙的方向与沉积层中的主要断层方向存在差异。后者对应于渐新世地堑展开的莱茵河走向N20°E，而花岗岩中大量裂隙方向为N160°E，与海西运动和阿尔卑斯山脉构造有关。从这方面看，裂隙区域的群集现象无法用断层来解释，但是Dezayes等［23］也没有得到其他解释。环测试这些方面对相关研究进行综述。

2.2 水力激发

水力激发可分为物理激发和化学激发这2种，前者通过水流的力学作用，使原已存在的密封裂隙重新张开，或者使靠近临界应力状态的岩体碎裂，而后者主要通过水流的化学作用对自然裂隙中的填充矿物进行酸性溶解。通常，在物理激发初期也会在注射水中添加一定浓度的重卤水，其目的是实现对井下深部高热区域进行激发［33］。Soultz试验场的多次水力激发试验表明，以卤水填充物的形式增加注射水的密度对水力激发非常重要，因为其强烈影响着热储层的初始状态以及生成(影响其向上、水平还是向下生长)［34］。另外，重卤水有助于水力激发产生的出水口沿着深度更均匀的分布，有效避免出水口过分集中于casing shoes附近。

目前，化学激发的方式是通过内套环往井头注入酸性液体，按照酸性液体的成分，化学激发可分为传统酸系统(conventional acid systems)，蛰合剂(chelatants)和缓速酸系统(retarded acid systems)这3种［35-36］。传统酸系统使用的酸有:(1)盐酸，用于溶解碳酸盐，如方解石和白云石等;(2)常规泥酸，用于溶解粘土和长石等矿质。蛰合剂使用氨三乙酸，通过形成阳离子复合物以减少阳离子的活性，从而增强矿石溶解作用。缓速酸系统使用有机黏土酸激发，用于一般酸质敏感的高温和高黏土含量条件，有机黏土酸的迟滞特性可以对深层热储进行激发。表1给出了GPK2，GPK3和GPK4井使用过的酸性激发情况［35］。

表1 各井使用过的化学激发情况Tab.1 Chemical stimulation methods for Soultz wells

从激发效果上看，虽然GPK2和GPK4井的总出水量提升了50%，但是Nami等［36］认为至少有一部分是套环水泥脱落的2个地方的物理激发贡献的，GPK3井就没有表现明显的改善，因此化学激发的作用还有待研究。Nami等［36］同时提出，化学激发的作用可以通过区域选择性的激发方式来提高，例如对断裂晶体结构中的裂隙区域接合处进行“聚焦”酸性激发。Portier等［37］认为化学激发和水力激发的协同可以减少水力激发对高水压的依赖，从而可以降低高压水力激发带来的地震危险。

2.3 微震监测

对微震云图演变的分析是研究水力激发条件下压裂事件的首要工具。为了对激发产生的微震进行监测和定位，Soultz地区布置了井下监测和地面监测这2种地震监测网。井下监测网由四轴加速计探头和三轴声波探头连接于高速数据采集仪和处理系统构成，探头分别置于4550，4601，EPS1，OPS4 和 GPK1这5口井底。地面监测网有35个站点，用于对较大的震动事件进行监测［33］。

2.3.1 数据处理

为了对水力激发的微震源进行定位，各学者采用了不同的数据处理方法对微震监测数据进行处理。通常采取的是由 Jones和 Stewart［38］提出来的collapsing方法，在假设震形结构比较简单的条件下得出其图形。联合震源定位(joint hypocentral determination，JHD)反演法［39］可用于提高震形结构图的精度，因此通常将collapsing方法作为约束条件与JHD反演法结合使用［40］。collapsing是一种倒转随机误差的统计优化方法，这种处理会导致相近的震源位置趋向分离，比较适用于移除震云结构中的模糊点，但对于小于误差椭圆的结构则无法处理［41］。为了获得表征微震活动云图结构的精微特征，一些学者使用了有别于基于统计论的collapsing方法而基于确定性的其他方法，例如双重/多重态分析法(doublet/multipletanalysis)［42］，群集分析法(clustering analysis)［43］以及多重－群集综合分析法(multipletclustering analysis)［44］等。这些方法都是通过减小P波和S波触发时间误差来提高震源定位精度的。此外，基于波形关联技术的相对迁移法［45］也在提升定位精度上具有较好的效果。

通过微震云图研究水力激发过程，不仅需要对震源准确定位，还需要将震源演变和流动参数同步对应，为此 Cuenot等［46］应用了三维局部断层扫描法(three-dimensional local tomography)实现。该方法通过同步转置可在三维各向异性介质中同时检索微震速度参量和震源定位参量:在定位计算的每次迭代中都考虑微震速度结构的演变。通过这种方法以及整个水力激发过程中微震事件的完整记录，Cuenot等［46］获得了震源演变和流动参数的同步比较分析。

2.3.2 微震分析

2000年6月30日，GPK2井进行了一次为期数日的水力激发试验，图4给出了微震事件和水力参数随时间的对应关系，图5给出了最终的微震云图［47］。本次激发按照不同的注水流量分为3个步骤，Step I和Step I、Step II的井头压力变化都有先迅速增大到峰值随后缓慢减小的共同特征，而Step III有所不同，井头压力持续增大直到注水关闭才缓慢减小，这种现象表明受激发区域是一个未与可以吸释压力的断层相连的封闭系统。激发过程中一共监测到了30 000 个微震事件，其中有14 000 个被 Dyer［48］进行了定位。Dorbath等［47］又对这些定位震源进行了筛选，确定了7 215个可信度较高的定位，并挑选其中M＞1(M为地震量级)的震源绘制了图5的微震云图。

图4 2000年GPK2井水力激发的微震事件和水力参数Fig.4 Seismic events and hydraulic parameters in 2000 stimulation of GPK2

图5 2000年GPK2井水力激发的最终微震云图Fig.5 Clouds of induced seismic events in 2000 stimulation of GPK2

图5显示了一个浓密、紧凑和均匀的震云，看不出明显的结构特征。为了研究激发过程中微震云图的演变情况，Dorbath等［47］将此次激发过程分为6个阶段:Phase I(30 L/s)，Phase II(40 L/s)，Phase III和Phase IV(分别对应Step III的前后两部分)，Phase V和Phase VI(分别对应注水关闭前、后)，并研究了这6个阶段的微震云图演变情况。研究发现，在激发过程中以及注水关闭之后的很长时间微震活动主要向北和向下发展，在Phase I和Phase II阶段微震云生长迅速，2天以后缓慢生长到最终结构。值得一提的是，在大约4 400 m深的地方(大概对应图3的Cluster III裂隙区域)发生了1次M=2.6的较大级微震事件。

2003年5月27日，GPK3井进行了一次为期数日的水力激发试验，图6给出了微震事件和水力参数随时间的对应关系，图7给出了最终的微震云图［47］。本次激发采用了双井注入方案，除了主要从GPK3井注入了33 000 m3的水量外，也从GPK2井短期注入了4 300 m3的水量。通过对比图4、6的微震活动情况可以发现，在2000年GPK2井注水关闭之前微震活动已经开始减少，而在2003年GPK3井注水关闭之后，微震活动似乎已成独立事件，不受水力参数的影响继续发生。Dorbath等［47］对这2次激发的b值(反映地壳强弱的参数)计算发现，GPK3井可能处于断层构造地带，一旦断层被激发，微震活动就只受断层本身机制作用而与注水计划无关。因此，从避免地震危害的角度考虑，未来EGS钻井的选址应当远离断层构造区域，而对这类区域的勘探是未来选址的重要任务。

图6 2003年GPK3井水力激发的微震事件和水力参数Fig.6 Seismic events and hydraulic parameters in 2003 stimulation of GPK3

图7 2003年GPK3井水力激发的最终微震云图Fig.7 Clouds of induced seismic events in 2003 stimulation of GPK3

图7的最终云图也显示出了北向发展的特征，与图5的浓密、紧凑和均匀特征不同的是，此次微震云图具有明显的发散特征。Baria等［33］对本次激发的水力参数和微震演变进行了详细研究，确定了4个裂隙结构(图中 F1、F2、F3、F4)。

F1位于4 660～4 710 m，具有偏西方向55°倾角，在注水最初阶段开始形成。F2位于4 870～4 920 m，具有偏西方向65°倾角，在F1之后随之形成，具有沿着井垂直向上生长的趋势。F3是一个较大的裂隙结构，呈向上和向下生长的趋势，也是偏西方向约65°倾角。值得注意的是，此次激发大部分较大的微震活动都发生在F3，而且其包含的4 400 m也与2000年GPK2井激发时M=2.6的较大微震事件发生点一致。F4在4 000 m处与F3相交，呈偏西方向40°倾角，向上东向生长与GPK2井在其3 870 m拐弯处相交。

2.4 循环测试

对注水试验和循环测试中的流量－压力曲线、流量记录、温度记录以及示踪浓度曲线的研究也极大地丰富了地下热储层的知识。1993年9月，Soultz EGS进行了几次重要的注水试验，图8给出了其井头流量－压力曲线［49］，从中可以看出注水过程中传输率的变化。岩石的渗透率一开始非常小(0.15 L/s阶段时为1.5×10－17m2)，当注水压力上升到5 MPa时开始迅速增大，表明了岩石中裂隙开始形成和生长，与微震监测的5 MPa监测到第1个微震事件的结果一致［50-51］。当流量上升到 6 L/s阶段，压力上升到7.6 MPa的稳定值，此时岩石的渗透率比初始值增大了20倍以上。此后当流量继续增大时，压力的增大逐渐减小，最终以9 MPa作为上限，表明裂隙的生长已经完成，热储层处于一个稳定的状态。

2.4.1 流量和温度研究

受井下裂隙分布结构影响，流量沿着深度变化的曲线会呈现出阶梯状，如图9所示，这些阶梯反映了流体流进或者流出岩石的位置，从而可以推断井下裂隙的分布情况［52］。以图9为例，最上面50 m包含5个离散的流体进出点，每个进出点都对应着独立的裂隙结构。再往下底层700 m只包含5个明显的阶梯，以100～150 m的距离间隔。基于这种分布规律，Evans［50-52］将井下这部分长度分为6个区域。但是单纯的流量记录只能看出主要的流体进、出点，而受到横截面(cross-sectional area，CSA)流体波动的噪声和阻力干扰因而无法对小流量的细微点进行探测。虽然CSA波动影响可以从理论上修正，但是实际上CSA修正能量的残差经常存在，这种问题在扰动剧烈的流体进出区域尤其严重。温度记录可以克服这些问题，从而可以精确判断流体进、出点，其原理是:由于更深层岩石的加热，井中上行的流体温度通常会高于周围岩石的温度，因而当有流体从岩石中流出时会造成井中流体温度的阶跃现象，井中流体下行时也会类似。

图8 Soultz EGS测试的井头流量随压力的变化曲线Fig.8 Chang curve of well head flow with pressure in Soultz EGS test

图9 GPK1井1993年注水试验的井下一段流量分布记录Fig.9 Selection of flow profiles from spinner logs during 1993 injection of GPK1

Evans［50］结合流量记录和温度记录对 GPK1井下一段长度35 m的裂隙结构进行了详细分析，如图10所示。从超声井像(UBI)分析判断，这一段包含2个先于水力激发的自然裂隙。从流量记录上看，这2个自然裂隙并非流体主要进出点，流体主要从高于他们的Z4-F1裂隙进出。除此之外，流量记录并没有显示此区域还有其他的流体通道，只能看到3 230～3 235 m间有流体进入，但是由于CSA波动影响而无法判断准确位置。这种问题在温度记录的辅助下得到了解决，大部分温度曲线都显示了2个阶梯，阶跃位置处于3 235 m处，正是其中一个自然裂隙所在，表明这一个自然裂隙有流体流通。

图10 GPK1井下一段长度35 m的超声井像、流量记录和温度记录Fig.10 Ultrasonic well image，flow and temperature record for a 35 m section of GPK1

2.4.2 示踪曲线研究

在历次循环试验中，Soultz EGS都在循环流体中掺入了一些化学物质作为示踪分析，这些示踪物质的种类 有 Na-benzoate，1，5-naphthalene disulfonate(nds)，2，7-nds，1，6-nds和2，6-nds以及荧光剂。对这些循环试验的示踪分析可以大致研究各井之间的连通情况，其中有重要意义的一次研究是2005年7月至12月的为期5个月的循环测试。向GPK3井以大约15 L/s的速率注入冷水，在开始阶段向水中加入150 kg 82.5%纯度的荧光剂，GPK2和GPK4井分别产出了11.9 L/s(80%)和3.1 L/s(20%)的热水，产出水的荧光剂的浓度曲线如图11所示［53］。荧光剂是在测试开始之后第8天掺入GPK3井的，4天之后即从GPK2井探测到，其浓度在第9至16天达到峰值(630～770 μg/s)并在接下来5个月期间逐渐减少。而GPK4井是在荧光剂注入之后第28天才开始探测到的，比GPK2井晚了24天，而且浓度也远低于GPK2井，在整个测试期间都没有达到峰值。从以上GPK2和GPK4井的产水量和荧光剂浓度可以看出，GPK3井与GPK2井之间有比较好的水力连通，而GPK3井与GPK4井之间的连通效果不佳。GPK2和GPK4井的荧光剂浓度变化趋势也引起了大家对热储层中水流路径的研究兴趣。

图11 2005年循环试验中GPK2井和GPK4井的荧光剂浓度随事件变化Fig.11 Fluorescein concentration change in fluid discharged from wells GPK2 and GPK4 in 2005 circulation tests

Sanjuan等［53］用离散传输模型分别对 GPK2和GPK4的荧光剂浓度数据进行了计算模拟，其中GPK2井获得了很好的吻合效果。由计算结果推断，GPK3和GPK2之间存在着2条连通路径，一条是短而快的循环通路Loop I，另一条是长而慢的循环通路Loop II，如图12所示。由于GPK3和GPK4的连通性较差，因此Sanjuan等［53］推断其间没有连通路径，而只在 Loop II有1条支线 Loop III通向 GPK4。Sanjuan等［53］关于 GPK3和 GPK2，GPK4井之间通路分析获得了后来众多学者的认同，他们都基于这种3通路结构对此次示踪试验进行了计算研究。Blumenthal等［54］采用简单二维数值模型，由于把流体路径简化在平行通道中，计算结果与流量曲线并不吻合。Genter等［55］在离散裂隙网络中应用粒子追踪方法构建计算模型，首次在同一个组合模型中对GPK2和GPK4井的浓度曲线进行了同时模拟，其中GPK2井的吻合较好。Kosack等［56］分别用全物理梯度贝叶斯转换法(full-physics gradient-based Bayesian inversion)、蒙特卡洛法(monte carlo approach)和总体卡尔曼滤波吸收法(ensemble-Kalman-filter assimilation)进行了计算，虽然全物理梯度贝叶斯转换法的计算结果与试验曲线有较好吻合，但是存在渗透率和孔隙率设置相反的不合理问题。Vogt等［57］也采了用蒙特卡洛法进行研究，与Kosack等［56］对渗透率均匀分配不同的是，Vogt等［57］用地球统计学对渗透率进行随机分配，使模型计算和实验数据获得了随机吻合。Radillaa等［58］采用分层多孔介质模型进行模拟，GPK2和GPK4井都取得了较好的吻合效果，为其通道构造给出了较好解释。

图12 GPK3井与GPK2井、GPK4井之间的连通情况推测Fig.12 Connectivity conjecture between GPK3 and GPK2，GPK4 wells

3 数值模拟

根据Soultz等野外试验场的研究，一个有效的热储层是由众多相互连通的裂隙带构成的复杂网络，如何构建热储层的数值模型并对热量传递过程(Thermal-T)、流体流动过程(Hydraulic-H)、力学过程(Mechanical-M)和化学过程(Chemical-C)进行THMC耦合计算是EGS数值模拟研究的重要问题。目前THMC多场耦合技术已经比较完善，一些学者都比较成功地研究了各个过程的耦合问题［59-62］。在热储层模型的构建方面仍然有待发展，将会继续成为EGS数值模拟研究的重点内容。此外，近来有学者提出了用超临界二氧化碳(supercritical CO2，SCCO2)代替水作为EGS换热工质的理念，并进行了初步的数值模拟研究。由于目前所有EGS试验场都没有用SCCO2工质的试验研究，因此数值模拟在EGS工质研究方面就比较重要。

3.1 模型研究

根据数值计算对热储层模型的处理方法，热储层计算模型大体上可以分为规则裂隙网络模型(图13［63］)、随机裂隙网络模型(图14［64］)和等价多孔介质模型等。规则裂隙网络模型是把热储层原本复杂的裂隙结构简化为规则排列的裂隙网格，包括单一裂隙(single fracture)和平行裂隙组(parallel fracture arrays)等，这种简化处理有助于用解析法对热储层中的热交换过程进行简单的理论分析。采用规则裂隙网络模型，DuTeau 和 Swenson［65］对 Fenton EGS 的岩石变形、流动和传热进行了耦合计算，发现岩石的冷却和变形对流体的传输阻抗有重要影响。Taron和Elsworth［62］在规则裂隙网络模型中增加裂隙的形变控制，研究了热储层中水力和化学激发过程，发现通过机械剪切产生的大尺度裂隙的激活是造成采热循环短路的重要因素。

图13 规则裂隙网络模型Fig.13 Regular fracture network model

一个实际的热储层包括裂隙方向、裂隙孔径、裂隙维度和裂隙密度等复杂信息，而规则裂隙网络模型无法将这些综合信息考虑在内，因此需要用随机网络模型来对热储层进行构建。随机裂隙网络模型是采用蒙特卡洛法确定热储层的平均特征，生成众多的结构各异但统计学上等效的裂隙网络。随机裂隙网络模型最先由Long［66］和一些学者在研究“migration of chemical species from waste disposals to bedrock”的问题中提出来的，Bruel和 Cacas［67］以及 Lanyon 和Batchelor［68］将其应用到EGS热储层模型构造问题。虽然随机裂隙网络模型比规则裂隙网络模型更接近热储层研究实际，但其存在由于自由变量过多而自由度太大的问题，给计算带来了沉重负担。为此，Watanabe 和 Takahashi［69］采用分形法来约束随机参数的变化。同时，为了限制计算任务，裂隙和岩体之间的换热也采用简化的传输函数或者一维扩散函数来近似处理，这样又给数值计算的有效性带来了问题［70］。

图14 随机裂隙网络模型Fig.14 Stochastic fracture network model

Bruel［71］在随机裂隙网络模型的基础上用确定性信息作为条件对计算进行约束，发展为确定裂隙网络模型。该模型仍然面临着由大量确定性信息带来的繁冗耗时的问题，而且Bruel［72］的计算结果也表明只有极少部分约束条件对结果有较大影响。Kolditz和Clauser［64］用实测数据结合一些假设对确定裂隙网络进行调整和简化，这些假设包括:(1)裂隙方向—流动传输限制在离散的裂隙系统中，这些裂隙的方向受地质构造应力的影响，所有流体路径都沿着平均走向和倾角;(2)裂隙分布—用大裂隙来简化裂隙束组，以平均流动简化分散流动;(3)裂隙间距—等距设置。这些调整和简化使得确定裂隙网络模型的计算量极大地减少，而且模型的有效性也在Rosemanowes EGS的循环采热算例中获得了验证。

多孔介质模型由于处理简便而广泛用于热储层模型的构建。一般的单孔隙模型在模拟从热储层吸热时都把岩石和流体之间假设为瞬态热平衡，计算结果一般会高估产热量。Shaik等［72］采用热非平衡模型，将岩石和流体之间换热系数考虑进计算中，强调了换热系数对热储层换热模拟的重要性。Yang和Yeh［73］以及Gelet等［74］采用了双孔隙率的方法来区分岩石和流道的物性，实现岩石和流体之间的热非平衡交换。为了更精确模拟热储层的各向异性效果，Kalinina等［75］提出了多孔隙率的多孔介质模型，将裂隙方向、间距、孔径和长度等重要裂隙参数都考虑进模型中。

如何在多孔介质模型中实现模拟岩石和流体之间的局部热交换一直是学者研究的问题，Jiang等［76］亦基于岩石骨架和流体之间的热非平衡交换假设、采用双能量方程模型实现了岩石和流体之间局部换热的模拟。

双能量方程模型如图15、16所示，该模型将EGS地下部分简化为3个主要区域:Region I为热储层，设置为具有一定孔隙率(ε)和渗透率(K)的多孔介质;Region II为无渗透的致密岩石，即设置为ε=0，K=0，隐含假设了流体损失为0;Region III为注入井和生产井，设置为ε=1和K=∞的开放通道。

工质水从注入井注入、采出井采出，与热储层发生热量交换，假设为单相和层流流动，且不考虑热储层形变和化学反应等作用。该模型所使用的控制方程如下。

图15 双能量方程计算模型地下区域划分Fig.15 Underground division in dual energy equations calculation model

图16 双能量方程计算模型尺寸Fig.16 Geometry of dual energy equations calculation model

质量守恒方程:

动量守恒方程:

岩石固体中的能量守恒方程:

孔隙流道中的能量守恒方程:

式中:ρ为密度;t为时间;u为速度矢量;P为压力;μ为动力粘度;k为渗透率;ε为孔隙率;g为重力加速度矢量;cp为定压比热容;T为温度;keff为有效渗透率;h为对流换热系数;a为换热面积;下标s和L分别表示固体和液体。式(2)中的第3项为达西阻力项。由于考虑岩石与流体之间的热非平衡过程，因此分别用式(3)、(4)来描述岩石固体中的导热和孔隙流道中的对流换热。其中 ±ha(Ts－Tl)用来描述岩石和流体之间的热量交换，是本模型区别于瞬态热平衡模型的关键设置。该模型已被用来模拟了多种工况下EGS的长期运行过程，分析了孔隙率、渗透率、循环流量、岩石－流体对流换热系数和热储裂隙比表面积等多个物理参数对采热的影响极其作用机制。图17示出了使用该模型模拟的不同渗透率热储中岩石温度的演化过程，图中上排的渗透率为1.0×10－12m－2，下排的为1.0×10－8m－2，时间由左向右递增。

图17 使用自主开发的模型模拟的不同渗透率下EGS运行过程中热储内岩石温度Fig.17 Rock temperature evolution with different permeability in self-developed model during EGS running process

3.2 工质研究

目前，EGS试验场的水力压裂和采热都是以水作为工质，而水资源的来源、水的损失以及水中矿物质的沉淀、岩石中矿物质的溶解等问题也是EGS研究和应用中需要考虑的。Brown［77］首次提出以SCCO2代替水作为EGS采工质的理念，指出了其三大优势:(1)注入井的冷 SCCO2和生产井的热SCCO2之间的显著密度差异能够提供强大的浮升力，从而减小工质循环所需的泵功;(2)SCCO2对矿物质具有极低的溶解性和传输能力，不会造成矿物质从热储层传输到地面设备的问题，因此对地面的管道、换热器和其他设备起到保护作用;(3)有利于温度高于374°C(水的临界温度)的热储层的开发，不会造成以水为工质时硅溶解导致的相关问题，从而提高热力效率。

Brown［77］用朗肯循环对以SCCO2作为采热工质的EGS电站进行了计算评估，发现SCCO2不仅具有和水一样好的热采能力，而且由于高比热容特性，在相同的注入井口压力、热储流动阻力和生产井口压力等条件下，SCCO2的密度和粘性之比(衡量热储层流动潜力的标准)是水的1.5倍。此外，从温室气体的处理方法来看，还可以把EGS工程与CO2封存工程相结合，以实现减少 CO2排放的目的［78］。Pruess［79］进一步全面地比较了SCCO2和水作为EGS采热工质的优点和劣势，如表2所示。

表2 二氧化碳和水作为EGS采热工质的比较Tab.2 Comparisons of CO2and water as work fluid of EGS

Pruess［80］首次对以 SCCO2作为工质的 EGS进行了数值模拟研究，评估了其系统的热学特性。他以一个近似Soultz EGS参数设置的理想热储层作为研究对象，构建了一个二维的1 km2区域的5井结构模型，如图18所示。这种结构模型由于高度的对称性而经常用于热储层特性的基础研究，给问题带来了极大的简化处理。整个计算范围只需要覆盖1/8的区域，其结果却可以对整个区域有效。该模型假设其中存在3组正交的裂隙组，相互间隔50 m，因此被裂隙包围的岩石骨架就是边长(略小于)50 m的正方形结构。岩石骨架视为无渗透能力，并进一步把网格细分为5个continua。虽然 Pruess［80］用该模型计算获得了一些结论，例如验证了工质CO2在采热方面优于工质水的能力，以及在维持循环方面少于工质水的耗功，但是该模型还存在2个较大的缺陷:(1)由于只耦合了换热和流动过程而未考虑化学作用，因此CO2和岩石矿物的相互作用以及CO2在化学特性方面的优势都没有研究;(2)CO2在减少循环耗功方面的优势来源于上下密度变化产生的浮升力作用，这是不考虑井深方向的二维平面模型无法进行研究的。为了考虑CO2密度变化产生的浮升力作用，Pruess［80］将他的二维模型扩展到三维，在原二维模型的基础上把热储层厚度分为6层，但总深度只有305 m。目前用SCCO2作为EGS工质的数值模拟研究仍然不够充分，未来的研究需要在三维模型和化学作用耦合上深入。

图18 二维5井结构模型和对称性处理Fig.18 Symmetrical simplification of 2-D five-well structure model

4 研究展望

通过野外试验场研究和数值模拟研究，EGS取得了一系列工程和理论成就。针对我国EGS研究起步晚、基础薄弱的研究现状，迫切需要开展如下3个方面的工作。

4.1 进行小型综合实验研究

建立野外试验场不仅需要投入大量的资金和技术力量，而且特定地区的地质特征和干热岩物性对实验结果影响巨大，使得不同地区EGS经验的类比和借鉴价值不大。相对于野外试验场的耗资耗时，实验室规模的小型综合实验系统具有操控性好、灵活性大、花费少等优势。比起野外试验场受制于特定地区的局限，小型综合实验系统可以方便地更换岩石材料进行研究，获得的成果具有更大地区范围的借鉴意义。而且由于实验成本小和灵活性高，可以为数值模拟提供更丰富的实验数据，促进数值模拟方法的完善和准确性的提高。可以预计，实验室规模的综合实验研究和数值模拟研究的紧密结合将会作为重要的研究方式出现在未来的EGS研究中。

4.2 建立多尺度数值模型

EGS地下热储层模型的尺寸通常在km量级甚至更大，而其内裂隙开度在mm量级甚至更小，裂隙分布十分不均匀且具有显著的各向异性。流体在狭小的裂隙中流动和换热，又在热储层尺度上构成整体循环，因此是涉及介观尺度和宏观尺度的多尺度过程。而且介观裂隙尺度的对流换热过程对宏观热储尺度的热开采过程具有决定性影响，如何解决介观尺度与宏观尺度之间的关系是开发热储层模型所要考虑的重要问题。在前文综述的各模型中，大多以规则网络法或等价多孔介质法来简化模拟热储层内复杂的裂隙网络结构，在模拟的真实性以及计算的准确度方面都存在不足。随机性和确定性裂隙网络模型虽然可以对裂隙微结构进行一定程度的真实构建，但也往往带来了网格数目及计算量过于庞大而无法计算的问题。解决以上问题的方案即是开发多尺度模型，即用介观模型对介观尺度过程进行研究，从其结果中提取更准确的参数应用到宏观模型中，既确保对热储层裂隙结构的准确构建，又使计算量控制在可接收范围。

4.3 野外试验场建设

建立野外试验场是研究EGS的最直接和最可信的方式。自从1974年美国在Fenton山从事干热岩试验研究以及随后一些发达国家相继建立野外试验场以来，EGS野外试验场已经有将近40年的研究历史，而我国目前尚未建立一座野外试验场。一个野外试验场的建设包含钻井探测、水力激发、微震监测、循环测试等几个阶段，完全可以借鉴石油天然气勘探开发的成熟技术，油气工程中的钻井技术、测井技术、水力压裂等技术可以用于EGS野外试验场的建设。我国的油气工程技术比较先进，而且众多学者已经在地热资源勘探和分析方面获得了众多成果，完全有条件开始建立野外试验场进行EGS研究。

5 结语

本文综述了国内外EGS研究的相关文献，介绍了国外野外试验场的钻井探测、水力激发、微震监测、循环测试等过程的研究成果和工程经验，也介绍了国内外数值模拟的模型构建、工质选择等问题的研究现状。我国的EGS研究还处于起步阶段，迫切地需要建立一个野外试验场来进行工程实际研究，同时还需要展开小型综合试验研究以及多尺度数值模型研究，提升我国EGS的发电技术和能力。

［1］国家能源局，财政部，国土资源部，住房和城乡建设部.关于促进地热能开发利用的指导意见［EB/OL］［2013－01－10］.http://zfxxgk.nea.gov.cn/auto87/201302/t20130207_1581.htm.

［2］张建英.增强型地热系统(EGS)资源开发利用研究［J］.中国能源，2011，33(1):29-37.

［3］汪集旸，胡圣标，庞忠和，等.中国大陆干热岩地热资源潜力评估［J］.科技导报，2012，30(32):25-31.

［4］蔺文静，刘志明，马峰，等.我国陆区干热岩资源潜力估算［J］.地球学报，2012，33(5):807-811.

［5］曾克棪，潘小平.中国地热发电开发现状与前景［J］.中外能源，2009，14(2):45-48.

［6］毛红举，马洪伟.增强型地热发电技术及广东省应用前景分析［J］.能源工程，2010(5):25-28.

［7］彭新明.我国干热岩地热资源钻采工艺浅议［J］.探矿工程，2009(S):167-169.

［8］赵阳升，万志军，康建荣.高温岩体地热开发导论［M］.北京:科学出版社，2004:23-69.

［9］李虞庚，蒋其垲，杨伍林.关于高温岩体地热能及其开发利用问题［J］.石油科技论坛，2007(1):28-40.

［10］杨吉龙，胡克.干热岩(HDR)资源研究与开发技术综述［J］.世界地质，2001，20(1):43-51.

［11］康玲，王时龙，李川.增强地热系统EGS的人工热储技术［J］.机械设计与制造，2008(9):141-143.

［12］Tester J W，Anderson B J，Batchelor A S，et al.The future of geothermal energy:Impact of enhanced geothermal systems(EGS)on the united states in the 21st century［R］.Massachusetts Institute of Technology，2006.

［13］Dezayes C，Genter A，Valley B.Structure of the low permeable naturally fractured geothermal reservoir at Soultz［J］.Geoscience，2010，342(7):517-530.

［14］Cornet F H，Jones R H.Field evidence on the orientation of forced water flow with respect to the regional principal stress directions［C］//1st North American Rock Mechanics Symposium.Texas，USA:Nelson P P，1994:61-69.

［15］Brudy M，Zoback M D.Drilling-induced tensile wallfractures:Implications for determination of in-situ stress orientation and magnitude［J］.International Journal of Rock Mechanics ＆ Mining Science，1999，36(2):191-215.

［16］Berard T，Cornet F.Evidence of thermally induced borehole elongation:A case study at Soultz，France［J］.International Journal of Rock Mechanics＆ Mining Science，2003，40(7-8):1121-1140.

［17］Valley B，Evans K F.Stress state at Soultz-sous-forets to 5 km depth from wellbore failure and hydraulic observations［C］//32nd Workshop on Geothermal Reservoir Engineering.Stanford，California，USA:Stanford University，2008.

［18］Genter A，Traineau H.Hydrothermally altered and fractured granite in an HDR reservoir in the EPS1 borehole，Alsace，France［C］//17th Workshop on Geothermal Reservoir Engineering.Stanford，California，USA:Stanford University，1992.

［19］Genter A， TraineauH.BoreholeEPS1(Alsace， France);preliminary geological results from granite core analyses for hot dry rock research［J］.Scientific drilling 1992(3):205-214.

［20］Genter A，Dezayes C.Fracture evaluation in GPK1 borehole by using FMI data［R］.BRGM Report 1993.

［21］Genter A， CuenotN， DezayesC， etal.How a better characterization of a deep crystalline reservoir can contribute to improve EGS performance at Soultz［C］//1st European Geothermal Review.Mainz，Germany:Geothermal Energy for Electric Power Production，2007:34-40.

［22］Valley B，Dezayes C，Genter A.Multiscale fracturing in the Soultz-sous-forets basement from borehole image analyses［C］//Proceedings EHDRA Scientific Conference.Soultz-sous-forets，France:EHDRA，2007:EP6-WP5.

［23］Dezayes C，Genter A，Valley B.Structure of the low permeable naturally fractured geothermal reservoir at Soultz［J］.Geoscience，2009，342(7):517-530.

［24］Dezayes C，Genter A，Valley B.Overview of the fracture network at different scales within the granite reservoir of the egs Soultz site(Alsace，France)［C］//Proceedings World Geothermal Congress.Bali，Indonesia:EHDRA，2010.

［25］Sausse J，Dezayes C，Genter A.Characterization of fracture connectivity and fluid flow pathways derived from geological interpretation and 3D modelling of the deep seated EGS reservoir of Soultz(France)［C］//33rd Workshop on Geothermal Reservoir Engineering.Stanford， California， USA: Stanford University，2008.

［26］Castera J，Dezayes C，Calcagno P.Large-scale 3D geological model around the Soultz site［C］//Proceedings of the EHDRA scientific conference.Soultz-sous-forets，France:EHDRA，2008.

［27］Sausse J，Dezayes C，Dorbath L，et al.3D model of fracture zones at Soultz-sous-forets based on geological data，image logs，induced microseismicity and vertical seismic profiles［J］.Comptes Rendus of Geoscience，2010，342(7):531-545.

［28］Dezayes C，Genter A.Large-scale fracture network based on Soultz borehole data［C］//Proceedings EHDRA Scientific Conference.Soultz-sous-forets，France:EHDRA，2008.

［29］Dezayes C，Genter A，Gentier S.Fracture network of the EGS Geothermal Reservoir at Soultz-sous-forêts(Rhine Graben，France)［C］//Geothermal Resources Council Transactions.California，USA:Palm Springs，2004.

［30］Genter A，Castaing C，Dezayes C，et al.Comparative analysis of direct(core)and indirect(borehole imaging tools)collection of fracture data in the hot dry rock Soultz reservoir(France)［J］.Journal of Geophysical Research，1997，102(B7):15419-15431.

［31］Hooijkaas G R，Genter A，Dezayes C.Deepseated geology of the granite intrusions at the Soultz EGS site based on data from 5 kmdeep boreholes［J］.Geothermics，2006，35(5 －6):484-506.

［32］Bruel D.Impact of induced thermal stress during circulation tests in an engineered fractured geothermal reservoir［J］.Oil＆ Gas Science and Technology，2002(5):459-470.

［33］Baria R， MicheletS， BaumgartnerJ， etal.Microseismic monitoring of the world's largest potential HDR reservoir［C］//29th Workshop on GeothermalReservoir Engineering.Stanford，California，USA:Stanford University，2004.

［34］Baumgartner J，Jung R，Gerard A，et al.The European HDR project at Soultz-sous-forets:stimulation of the second deep well and first circulation experiments［C］//21st Workshop on Geothermal Reservoir Engineering.Stanford，California，USA:Stanford University，1996.

［35］Portier S，Andre L，Vuataz F.Review of chemical stimulation techniques in oil industry and applications to geothermal systems［R］.Technical Report，ENGINE work package 4:Drilling，stimulation and reservoir assessment， Neuchatel，Switzerland，2007.

［36］Nami P， SchellschmidtR， SchindlerM， etal.Chemical stimulation operationsforreservoirdevelopmentofthe deep crystalline HDR/EGS system at Soultz-sous-forets(France)［C］//33rd Workshop on Geothermal Reservoir Engineering.Stanford，California，USA:Stanford University，2008.

［37］Portier S， VuatazF， Nami P， et al.Chemical stimulation techniques for geothermal wells:experiments on the 3-well EGS system at Soultz-sous-forets［J］.Geothermics，2009，38(4):349-359.

［38］Jones R H，Stewart R C.A method for determining significant structures in a cloud of earthquakes［J］.Journal of Geophysical Research，1997，102(B4):8245-8254.

［39］Pujol J.Joint hypocentral location in media with lateral velocity variations and interpretation of the station corrections［J］.Physics of Earth Planetary Interiors，1992，75(1-3):7-24.

［40］Fehler M C，Phillips W S，House L S，et al.Improved relative locations of clustered earthquakes using constrained multiple event location［R］.Bulletin of the Seismological Society of America，2000，90:775-780.

［41］Evans K F，Moriya H，Niitsuma H，et al.Microseismicity and permeability enhancementof hydrogeologic structures during massive fluid injections into granite at 3 km depth at the Soultz HDR site［J］.Geophysical Journal International，2005，160(1):388-412.

［42］Poupinet G，Ellsworth W L，Frechet J.Monitoring velocity variations in crust using earthquake doublets:An application to the Calaveras fault［J］.Journal of Geophysical Research，1984，89(B7):5719-5731.

［43］Phillips S，House L S，Fehler M C.Detailed joint structure in a geothermal reservoir from studies of induced microseismic clusters［J］.JournalofGeophysicalResearch，1997，102(B6):11745-11763.

［44］Moriya H，Niitsuma H，Baria R.Multiplet-clustering analysis for estimation of fine detail structures in seismic clouds at Soultz field，France［R］.Bulletin of the Seismological Society of America，2003，93:1606-1620.

［45］Moriya H，Nakazato K，Niitsuma H，et al.Detailed fracture system ofthe Soultz-sous-foretsHDR field evaluated using microseismic multiplet analysis［J］.Pure and Applied Geophysics，2002，159(1－3):517-541.

［46］Cuenot N，Dorbath C，Dorbath L.Analysis of the microseismicity induced by fluid injections at the EGS site of Soultz-sous-forets(Alsace，France):Implications for the characterization of the geothermal reservoir properties［J］.Pure and Applied Geophysics，2008，165(5):797-828.

［47］Dorbath L，Cuenot N，Genter A，et al.Seismic response of the fractured and faulted granite of Soultz-sous-forets(France)to 5 km deep massive water injections［J］.Geophysical Journal International，2009，177(2):1-23.

［48］Dyer B C.Soultz GPK2 stimulation June/July 2000［R］.Seismic Monitoring Field Report，2000.

［49］Jung R， Willis-Richards J， Nicholls J， et al.Evaluation of hydraulic tests at Soultz-sous-forets，European hot dry rock site［C］//World Geothermal Congress.Florence，Italy:Barbier E，1995:2671-2676.

［50］Evans K F，Genter A，Sausse J.Permeability creation and damage due to massive fluid injections into granite at 3.5 km at Soultz，Part 1:Borehole observations［J］.Journal of Geophysical Research，2005，110(B4):1-19.

［51］Evans K F.Permeability creation and damage due to massive fluid injections into granite at 3.5 km at Soultz，Part 2:Critical stress and fracture strength［J］.Journal of Geophysical Research，2005，110(B4):1-14.

［52］Evans K F.The effect of the 1993 stimulations of well GPK1 at Soultz on the surrounding rock mass:evidence for the existence of a connected network of permeable fractures［C］//Proceedings of the World Geothermal Congress.Kyushu-Tohoku，Japan:Iglesias E，2000:3695-3700.

［53］Sanjuan B，Pinault J，Rose P，et al.Tracer testing of the geothermal heat exchanger at Soultz-sous-forêts(France)between 2000 and 2005［J］.Geothermics，2006，35(5 －6):622-653.

［54］Blumenthal M，Kühn M，Pape H，et al.Numerical modeling of hydraulic and thermal processes in the deep heat reservoir of the EGS test site Soultz-sous-forêts，France［C］//Proceedings of the Euro-Conference of Rock Physics and Geomechanics on Natural Hazards:Thermo-Hydro-Mechanical Processes in Rocks.Erice，Italy:EGU，2007.

［55］Gentier S，Rachez X，Ngoc T T，et al.3D flow modeling of the medium-term circulation test performed in the deep geothermal site of Soultz-sous-forets(France)［C］//Proceedings ofWorld Geothermal Congress.Bali，Indonesia:BRGM，2010.

［56］Kosack C，Vogt C，Rath V，et al.Stochastic estimates of the permeability field of the Soultz-sous-forêts geothermal reservoircomparison of bayesian inversion，MC geostatistics，and EnKF assimilation［C］//Proceedings of the European Geosciences Union General Assemble.Vienna，Austria:EGU，2010.

［57］Vogt C，Kosack C，Marquart G.Stochastic inversion of the tracer experiment ofthe enhanced geothermal system demonstration reservoir in Soultz-sous-forêts:Revealing pathways and estimating permeability distribution［J］.Geothermics，2012，42:1-12.

［58］Radillaa G，Sausse J，Sanjuan B，et al.Interpreting tracer tests in the enhanced geothermal system(EGS)of Soultz-sous-forêts using the equivalent stratified medium approach［J］.Geothermics，2012，44:43-51.

［59］Kohl T，Megel T.Coupled hydro-mechanical modelling of the GPK3 reservoir stimulation at the European EGS site Soultz-sousforets［C］//30th Workshop on Geothermal Reservoir Engineering.Stanford，California，USA:Stanford University，2005.

［60］Kohl T，Megel T.Predictive modeling of reservoir response to hydraulic stimulations at the European EGS site Soultz-sous-forets［J］.International Journal of Rock Mechanics＆ Mining Sciences，2007，44(8):1118-1131.

［61］Taron J，Elsworth D，Min K B.Numerical simulation of thermalhydrologic-mechanical-chemical processes in deformable，fractured porous media［J］.International Journal of Rock Mechanics ＆Mining Sciences，2009，46(5):842-854.

［62］Taron J，Elsworth D.Coupled mechanical and chemical processes in engineered geothermal reservoirs with dynamic permeability［J］.International Journal of Rock Mechanics＆Mining Sciences，2010，47(8):1339-1348.

［63］Hicks T W，Pines R J，Willis-Richards J，et al.A hydro-thermomechanicalnumericalmodelforHDR geothermalreservoir evaluation［J］.International Journal of Rock Mechanics ＆ Mining Sciences，1996，33(5):499-511.

［64］Kolditz O，Clauser C.Numerical simulation of flow and heat transfer in fractured crystalline rocks:application to the hot dry rock site in Rosemanowes(U.K.)［J］.Geothermics，1998，27:1-23.

［65］DuTeau R，Swenson D V.Modeling flow in a jointed geothermal reservoir［J］.Geothermal Resources Council Transactions，1994，18(6):465-472.

［66］Long J C S.Investigation ofequivalentporous medium permeability in networks of discontinuous fractures［D］.Berkeley，CA:University of California，1983.

［67］BruelD， Cacas M C.Numericalmodelling technique:contribution totheSoultzHDR project［M］.New York:Geothermal Energy in Europe—the Soultz hot dry rock Project，Gordon ＆ Breach，1992:57-74.

［68］Lanyon G W，Batchelor A S.Modelling for the UK hot dry rock programme using the fracture network method:Summary report［R］.Geoscience Report CC59/R5/078，1992.

［69］Watanabe K，Takahashi H.A model analysis of the long-term performance of hot dry rock geothermal energy extractions systems［C］//Proceedings of the JSME-ASME International Conference on Power Engineering.Tokyo，Japan:JSME-ASME，1993:453-458.

［70］Kolditz O.Modelling flow and heat transfer in fractured rocks:dimensional effect of matrix heat diffusion［J］.Geothermics，1995，24:421-437.

［71］Bruel D.Heat extraction modelling from forced fluid flow through stimulated fractured rock masses:application to the Rosemanowes hot dry rock reservoir［J］.Geothermics，1995，24:361-374.

［72］Shaik A R，Rahman S S，Tran N H，et al.Numerical simulation offluid-rock coupling heat transfer in naturally fractured geothermal system［J］.Applied Thermal Engineering，2011，31(10):1600-1606.

［73］Yang Y，Yeh H.Modeling heat extraction from hot dry rock in a multi-well system［J］.Applied Thermal Engineering，2009，29(8-9):1676-1681.

［74］Gelet R，Loret B，Khalili N.A thermo-hydro-mechanical coupled model in local thermal non-equilibrium for fractured HDR reservoir with double porosity［J］.Journal of Geophysical Research，2012，117(7):205-228.

［75］Kalinina E，McKenna S，Hadgu T，et al.Analysis of the effects of heterogeneity on heat extraction in an EGS represented with the continuum fracture model［C］//37th Workshop on Geothermal Reservoir Engineering.Stanford， California， USA:Stanford University，2012.

［76］Jiang F，Luo L，Chen J.A novel three-dimensional transient model for subsurface heat exchange in enhanced geothermal systems［J］.International Communications in Heat and Mass Transfer，2013，41(1):57-62.

［77］Brown D.A hot dry rock geothermal energy concept utilizing supercritical CO2instead of water［C］//25th Workshop on Geothermal Reservoir Engineering.Stanford，California，USA:Stanford University，2000.

［78］Pruess K.On production behavior of enhanced geothermal systems with CO2as working fluid［J］.Applied Thermal Engineering，2008，49(6):1446-1454.

［79］Fouillac C，Sanjuan B，Gentier S，et al.Could sequestration of CO2be combined with the development of enhanced geothermal systems?［C］//The third annual conference on carbon capture and sequestration，2004.

［80］Pruess K.Enhanced geothermal systems(EGS)using CO2as working fluid-a novel approach for generating renewable energy with simultaneous sequestration of carbon［J］.Geothermics，2006，35:351-367.