增強型地熱系統研究綜述

郭劍，陳繼良，曹文炅，蔣方明

(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.