雷宏武，金光榮，李佳琦，石 巖，馮 波

吉林大學(xué)地下水資源與環(huán)境教育部重點實驗室，長春 130021

0 引言

干熱巖（hot dry rock，HDR）作為地?zé)豳Y源的一種，由于分布廣泛、清潔和熱儲量巨大，被認(rèn)為是21世紀(jì)最有潛力的資源［1］。針對內(nèi)部不存在流體或僅有少量地下流體且滲透率極低的高溫巖體（干熱巖），增強型地?zé)嵯到y(tǒng)（enhanced geothermal system，EGS）或工程型地?zé)嵯到y(tǒng)（engineered geothermal system）作為一種新型的地?zé)岚l(fā)電技術(shù)，并不需要天然熱液資源的存在，可通過水力壓裂的方式從干熱巖中提取地?zé)豳Y源［2］。在地?zé)豳Y源的提取過程中，溫度場、壓力場、力學(xué)場和化學(xué)場互相影響，并隨著時空發(fā)生復(fù)雜的變化。采用傳統(tǒng)的方法無法全面了解這種復(fù)雜耦合過程的互相影響，而數(shù)值模擬作為地質(zhì)資源和環(huán)境研究的新方法［3］，在耦合過程的研究中發(fā)揮著越來越大的作用。Yano和Ishido［4］采用數(shù)值模擬方法研究了深部地?zé)嵩诔R界條件下的生產(chǎn)行為，表明壓力在超臨界區(qū)域變化復(fù)雜，其主要受流體動力黏滯性和壓縮能力的非線性特征控制。FRACTure和CHEMTOUGH被耦合形成FRACHEM，用于研究了Soultz深部地?zé)岬亩鄨鲴詈蠁栴}［5－7］。Xu等［8］采 用 TOUGHREACT研究了Soultz深部地?zé)嵯到y(tǒng)注入不同成分水的地球化學(xué)演化特征，以及其對注入能力的影響。Pruess［9－10］采用 TOUGH2對比分析了 CO2和水作為載熱介質(zhì)時的地?zé)嵯到y(tǒng)中的壓力和溫度變化特征，以及熱提取效率差別。

筆者以松遼盆地徐家圍子為研究對象，采用數(shù)值模擬方法研究了其深部地?zé)衢_發(fā)過程中溫度和壓力的變化特征，并分析了模型存在的各種不確定性，以期為下一步工程實施提供理論依據(jù)。

1 松遼盆地地質(zhì)概況

松遼盆地位于中國東北部，總面積約26萬km2，是一個大中型中新生代陸相沉積的、呈北北東向展布的菱形盆地。盆地主要有6個一級構(gòu)造單元［11］（圖1）：中央凹陷區(qū)、東北隆起區(qū)、東南隆起區(qū)、西南隆起區(qū)、西部斜坡區(qū)、北部傾沒區(qū)。盆地主要沉積了中、新生代的沉積巖系，最大沉積厚度大于10 km，地層從下而上依次為上侏羅統(tǒng)、白堊系、第三系和第四系。盆地深層是指泉二段以下地層，自上之下由泉頭組、登婁庫組、營城組、沙河子組和火石嶺組組成［12］。

松遼盆地現(xiàn)今平均地溫梯度是3.7℃/100m，是中國大型盆地中地溫梯度比較高的地區(qū)之一［13］。徐家圍子斷陷位于中央凹陷區(qū)東北部，主要由徐西斷裂、宋西斷裂和徐東斷裂分割成2個較大的斷裂區(qū)（圖2），是松遼盆地深部重要的油氣儲層，并且其靠近相對較高的地溫梯度帶，是優(yōu)良的深部地?zé)衢_采靶區(qū)之一。

宋西斷裂和徐東斷裂構(gòu)成的斷裂區(qū)沉積厚度較大，約4.5km，根據(jù) MIT報告［1］對增強地?zé)嵯到y(tǒng)的定義，適用于深部地?zé)衢_采的深度為3.0～10.0 km，因此，本區(qū)營城組中下部、沙河子組和火石嶺組將是深部地?zé)衢_采的備選層位。營城組為火山巖－陸源碎屑巖沉積建造，下部主要由安山玄武巖、火山角礫巖、砂巖及砂礫巖組成，上部為酸性火山巖、火山碎屑巖及砂巖和粉砂巖；沙河子組以泥巖為主，夾有砂巖、粉砂巖及少量凝灰?guī)r；火石嶺組主要以凝灰質(zhì)角礫巖、凝灰?guī)r、安山巖、玄武巖及凝灰質(zhì)礫巖為主。

圖2 徐家圍子斷陷地質(zhì)剖面概化圖（據(jù)文獻［12］修改）Fig.2 Geologic profile across the Xujiaweizi depression in E－W direction（modified after reference［12］）

2 模型建立

2.1 模擬工具

用于地?zé)衢_采的數(shù)值模擬軟件可分為2類：干熱巖和水熱型地?zé)嵯到y(tǒng)模擬器。干熱巖中典型的模擬 軟 件 有 FRACTure［14］、 GEOTH3D［15］、FRACSIM－3D［16］和 Geocrack2D［17］等，水熱型中常用的有 TOUGH2［18］、STOMP［19］、和 FEHM［20］等。TOUGH系列由于其能夠刻畫多個過程和可信的模擬結(jié)果［21］，被廣泛用于與地?zé)嵊嘘P(guān)的研究中［8，22－23］，其靈活的網(wǎng)格剖分方法和多重連續(xù)介質(zhì)模型［24］（multiple interacting continua，MINC）使得其不僅適用于水熱型地?zé)嵯到y(tǒng)，而且也能夠用于采用水作為載熱流體的孔隙和裂隙共同存在的干熱巖地?zé)嵯到y(tǒng)中。因此，本次研究擬采用TOUGH系列軟件進行。

TOUGH2為TOUGH系列的主要成員，其基于組分質(zhì)量和能量守恒求解多相流動過程中壓力、溫度、飽和度和組分質(zhì)量分?jǐn)?shù)的時空分布，采用的空間離散方法是積分有限差，時間離散是隱式差分，非線性方程求解是Newton－Raphson迭代［18］。TOUGH2采用模塊化設(shè)計，EOS1模塊能夠刻畫純水、水蒸氣及兩相共存系統(tǒng)，溫度最高可達647.3 K，其兩相飽和線、密度、黏度、熱焓等根據(jù)國際公式化委員會提供的模型計算。因此，TOUGH2/EOS1可用于松遼盆地深部地?zé)衢_采的模擬工作，其主要特征見文獻［18］。生產(chǎn)井和注入井的刻畫采用生產(chǎn)指數(shù)模型［25］。

2.2 模型概化

深部熱儲層孔隙度和滲透率一般都較低，往往需要進行人工壓裂形成新的裂隙來增強系統(tǒng)的流動能力。裂隙系統(tǒng)提供了流體快速流動的通道，保證了熱量可以循環(huán)不斷的提取。徐家圍子斷陷區(qū)深層為沉積巖和火山巖混雜，具有一定的孔隙度和滲透率，可以認(rèn)為是連續(xù)的孔隙系統(tǒng)。因此，通過人工改造以后的徐家圍子深部地?zé)嵯到y(tǒng)是一個孔隙－裂隙復(fù)合的復(fù)雜系統(tǒng)。

本次研究將徐家圍子深部地?zé)嵯到y(tǒng)在垂向上概化為孔隙和裂隙2個互相關(guān)聯(lián)的統(tǒng)一系統(tǒng)。假設(shè)壓裂范圍垂向上為50m，采用MINC刻畫。熱提取過程中由于上下孔隙系統(tǒng)的存在，注入井附近流體會向孔隙系統(tǒng)漏失，生產(chǎn)井附近流體會向裂隙系統(tǒng)流入，同時孔隙系統(tǒng)的熱會向裂隙系統(tǒng)進行傳導(dǎo)。為了刻畫這種影響，裂隙系統(tǒng)上下各取500m的孔隙介質(zhì)。平面上為了簡化邊界的影響，采用“五點法”布井方式［9－10］（圖3a）。由于對稱性，模型只需取1/8進行計算。平面上模型網(wǎng)格剖分間距為70.7m×70.7m，共剖分36個網(wǎng)格；垂向上網(wǎng)格間距為50.0 m，共21層。裂隙系統(tǒng)網(wǎng)格為180個，孔隙系統(tǒng)為720個，模型共900個網(wǎng)格（圖3b）。

徐家圍子斷陷深部可供地?zé)衢_采的層位有3層，為了分析開采過程中不同層位熱和水動力的異同，筆者建立了3個典型位置的模型（埋深為3 250、3 750、4 100m，分別以 Case1、Case2和 Case3代表），其對應(yīng)的中間壓裂層的溫度分別為131、147和160℃，初始壓力符合靜水壓力，整個模型的初始溫度按照T＝111＋0.035（D－2 500）（D＞2 500m）給定。注入井和生產(chǎn)井均采用定壓力方式，分別高于和低于初始壓力2MPa，模擬開采時間為36a。模型具體參數(shù)見表1。

圖3 “五點法”地?zé)衢_采概念模型和網(wǎng)格剖分Fig.3 Schematic representation of the five－spot well pattern geothermal development with computational grid

表1 “五點法”開采熱儲層模型基本物理參數(shù)Table 1 Basic property of geothermal storage formation in five－spot well pattern

3 模擬結(jié)果及分析

3.1 壓力場時空變化特征

從壓力時空化特征（圖4）可以看到：裂隙系統(tǒng)的壓力相對初始壓力大部分是降低的，生產(chǎn)井附近開采引起的壓力降低較大，而注入井由于持續(xù)注入壓力降低較?。徽麄€裂隙系統(tǒng)的壓力分布明顯受生產(chǎn)井的控制，但隨著時間的增加，注入的影響逐漸增大，整個系統(tǒng)壓力逐漸上升。這種變化特征主要原因是水的移動性（ρ/μ，其中：ρ和μ分別為流體密度和黏度）在注入和生產(chǎn)區(qū)域的不同，注入?yún)^(qū)域由于冷水的注入使得水的溫度降低，導(dǎo)致其移動性比生產(chǎn)區(qū)域的降低了約80%，從而使注入?yún)^(qū)域的壓力梯度明顯大于生產(chǎn)區(qū)域。隨著熱能逐漸被開采，生產(chǎn)井附近溫度降低，水的移動性的差別減小，注入的影響增大。

從圖4不同埋深位置的結(jié)果可以看到，埋深越大，整個系統(tǒng)的壓力降低程度越大，其主要原因是埋深大，溫度高，水的移動性增強，定壓產(chǎn)出條件下的水的流出質(zhì)量大，從而使得壓力降低大。

3.2 溫度場時空變化特征

注入的冷水在由注入井到生產(chǎn)井流動的過程中，與其接觸的介質(zhì)進行對流傳導(dǎo)換熱，由于裂隙介質(zhì)體積比例遠遠小于孔隙介質(zhì)，因此，換熱過程主要發(fā)生在流體和孔隙介質(zhì)之間?？紫督橘|(zhì)的低孔低滲特征又使換熱方式以傳導(dǎo)為主。

溫度變化由注入井向生產(chǎn)井是逐漸減小的；隨著時間的增加，由生產(chǎn)井到注入井的地?zé)嶂饾u提取完畢（圖5）。由不同埋深的模擬結(jié)果（圖5）對比可以看到：隨埋深增大，溫度變化增大，主要原因是埋深大的熱儲層初始溫度較高；另外其溫度變化影響范圍要大于埋深小的，主要由于溫度高加快了水的流動。

從圖6中孔隙基質(zhì)溫度的變化特征可以看到：與流體溫度變化特征相似，由生產(chǎn)井到注入井溫度逐漸降低；并且孔隙基質(zhì)中溫度分布存在梯度，隨著與裂隙距離的增加，溫度降幅減小，反映了孔隙基質(zhì)中由內(nèi)向外的傳熱過程。孔隙基質(zhì)中的這種溫度梯度在整個溫度變化的前緣較大，而在后緣逐漸變小，反映了后緣部分孔隙基質(zhì)中的熱逐漸被提取完畢。埋深大的溫度變化大，同樣是由于孔隙基質(zhì)本身的初始溫度較高，加快了水的流動速度。裂隙系統(tǒng)上下未壓裂處理的孔隙基質(zhì)的溫度變化不大，主要變化集中在壓裂區(qū)域上下附近，并且上下變化量很相近，說明未壓裂和壓裂區(qū)域之間的熱交換以傳導(dǎo)為主，對流傳熱可以忽略。

圖4 不同埋深位置地?zé)衢_采過程壓力變化時空分布Fig.4 Spatial and temporal distribution of pressure at different depth levels

圖5 不同埋深位置地?zé)衢_采裂隙介質(zhì)中溫度變化時空分布Fig.5 Spatial and temporal distribution of temperature in fractured media at different depth levels

圖7顯示了生產(chǎn)井和注入井連線（最短流動路徑）的溫度和壓力變化?？梢钥吹剑鹤⑷刖畢^(qū)域壓力上升范圍有限，大約50m，這主要是受注入井低溫的影響；注入井附近壓力梯度明顯大于生產(chǎn)井區(qū)域；不同埋深情況的壓力變化并不是很大，但溫度變化較大。

從圖8可以看出：生產(chǎn)井附近壓力由于生產(chǎn)井壓力降低而迅速降低并穩(wěn)定，在開采后期整個系統(tǒng)的溫度降低，溫度引起的流動阻力增大，從而壓力略有上升；埋深大，初始溫度和壓力均較高，從而在開發(fā)過程中其也相對較高。

3.3 熱提取速率

從圖9可以看到：對于埋深3 250m（Case1），由于開采過程中溫度逐漸降低，注入質(zhì)量速率和生產(chǎn)質(zhì)量速率由開始的70kg/s左右逐漸減小到53kg/s左右，并且注入質(zhì)量速率稍大于生產(chǎn)質(zhì)量速率；埋深越大，注入質(zhì)量速率和生產(chǎn)質(zhì)量速率均增大，相對埋深3 250m而言，埋深3 750m（Case2）和4 100m（Case3）的速率分別增大了2%～4%和4%～6%。埋深3 250m的熱提取速率由開始的40MW左右減小到15MW左右（圖10）。與質(zhì)量速率變化一樣，熱提取速率也隨著埋深增大而增加，埋深3 750 m和4 100m的熱提取速率分別增加了11%～17%和19%～29%，明顯大于質(zhì)量速率的增幅，主要原因是埋深大的流體溫度高，本身攜帶的能量相對較大。

因此，從熱提取速率的角度來說，松遼盆地徐家圍子3個埋深水平中埋深4 100m水平明顯要優(yōu)于其他埋深水平，同時也說明了埋深越大，溫度越高，熱提取速率也高。但深度增加，也必然增加了鉆井成本。

4 參數(shù)不確定分析

地質(zhì)條件的復(fù)雜性使得參數(shù)存在很多的不確定性，從而導(dǎo)致模擬結(jié)果的不確定性。為了評價這種不確定性，筆者進行了孔隙基質(zhì)滲透率和孔隙度、裂隙介質(zhì)滲透率和裂隙度、巖石導(dǎo)熱系數(shù)、井徑、注入/生產(chǎn)壓力、注入溫度等因素的不確定性分析，具體模擬設(shè)置見表2。

圖6 不同埋深位置地?zé)衢_采注入井－生產(chǎn)井之間孔隙介質(zhì)中溫度變化時空分布Fig.6 Spatial and temporal distribution of temperature in porous media at different depth levels

圖7 注入井－生產(chǎn)井之間壓力和溫度變化量分布（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 生產(chǎn)井所在網(wǎng)格壓力和溫度時間變化特征Fig.8 Evolution of pressure and temperature of production well

圖9 注入井和生產(chǎn)井質(zhì)量速率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 孔隙介質(zhì)滲透率和孔隙度

從圖11和圖12可以看出：孔隙基質(zhì)的滲透率和孔隙度對壓力和溫度的影響很小，整體上在靠近生產(chǎn)井區(qū)域壓力隨著孔隙基質(zhì)滲透率和孔隙度的增大稍微有所降低，溫度隨孔隙基質(zhì)滲透率和孔隙度的增大稍微有所升高；質(zhì)量和熱提取速率整體上是隨著孔隙基質(zhì)滲透率和孔隙度的增大而升高的。

4.2 裂隙介質(zhì)滲透率和裂隙度

從圖13可以看出：裂隙介質(zhì)滲透率對溫度、壓力、質(zhì)量和熱提取速率的影響很大；裂隙介質(zhì)滲透率逐漸增大時，注入井逐漸占主導(dǎo)地位，壓力增大，溫度降低；質(zhì)量和熱提取速率隨著裂隙介質(zhì)滲透率的增大而增大，當(dāng)裂隙滲透率增加到10－12數(shù)量級時，熱提取速率變化劇烈，由開始的160MW到開采15a時的10MW，熱儲存中的熱幾乎被開采完畢。因此，裂隙介質(zhì)滲透率直接決定了壓力和溫度的時空分布特征，以及質(zhì)量和熱提取速率。

表2 不同條件的數(shù)值模擬設(shè)置Table 2 Numerical simulation runs for different conditions

從圖14可以看出：裂隙介質(zhì)裂隙度對溫度、壓力、質(zhì)量和熱提取速率沒有影響。這是因為裂隙介質(zhì)所占比例很小，裂隙度的變化帶來的影響可以忽略，但是裂隙度直接影響流體的實際流速，從而可以改變整個系統(tǒng)的循環(huán)周期。

4.3 熱傳導(dǎo)系數(shù)

從圖15可以看出：壓力幾乎不受巖石的熱傳導(dǎo)系數(shù)的影響；溫度降幅隨著熱傳導(dǎo)系數(shù)增加而減小；而質(zhì)量和熱提取效率相應(yīng)地增加。這是因為熱傳導(dǎo)系數(shù)增加使得基質(zhì)中的熱更快地傳導(dǎo)到裂隙介質(zhì)中來，冷水注入引起的溫度降低減小，加快了流體的流動過程，從而提高了質(zhì)量和熱提取速率。但因巖石的熱傳導(dǎo)系數(shù)間差別較小，這種影響并不明顯。

圖11 孔隙介質(zhì)滲透率對溫度、壓力和質(zhì)量及熱提取效率的影響Fig.11 Effects of permeability of porous media on temperature，pressure，and mass and heat extraction rate

4.4 井徑

從圖16可以看出：隨著井徑的增加，井與地層的接觸面積增大，使得注入速率增大，導(dǎo)致壓力降低程度減小，溫度在靠近生產(chǎn)井區(qū)域降低程度增大。在開采過程中，隨著井徑的增大，質(zhì)量速率增大，而熱提取速率在開采早期增大，開采后期這種差別逐漸減小。這主要是因為井徑大的開采后期時系統(tǒng)溫度降低快，故單位質(zhì)量的流體攜帶的熱量變小。

4.5 注入井/生產(chǎn)井壓力和溫度

圖12 孔隙基質(zhì)孔隙度對溫度、壓力和熱提取效率的影響Fig.12 Effects of porosity of porous media on temperature，pressure，and mass and heat extraction rate

從圖17可以看出：注入井/生產(chǎn)井壓力大，在注入井區(qū)域（x＝0.0m）壓力增大得多，在生產(chǎn)井區(qū)域（x＝700.0m）壓力降低得多，生產(chǎn)井區(qū)域的變化幅度明顯大于注入?yún)^(qū)域，生產(chǎn)井控制整個系統(tǒng)的壓力，同時開采質(zhì)量和熱提取速率也高，溫度降低程度大；注入/生產(chǎn)壓力大的，隨著熱被迅速提取，溫度降低，其熱提取速率迅速變小。

從圖18a可以看出：注入溫度的增加可以增強注入?yún)^(qū)域水的流動性，提高注入井的影響范圍，從而使得整個系統(tǒng)的壓力降低較小，大大提高質(zhì)量流動速率；但是注入的能量明顯增加，最終的能量凈提取效率并不一定增加。圖18b顯示隨著注入溫度的增加，質(zhì)量速率增加，但熱提取速率整體減小。溫度分布在注入井附近注入溫度越大，降低越小，而在生產(chǎn)井附近由于注入溫度高的質(zhì)量速率大，從而其溫度降低得快。

圖13 裂隙介質(zhì)滲透率對溫度、壓力和熱提取效率的影響Fig.13 Effects of permeability of fractured media on temperature，pressure，and mass and heat extraction rate

4.6 基質(zhì)

從圖19可以看出：不考慮壓裂區(qū)域外孔隙介質(zhì)的影響，溫度降低明顯增大，導(dǎo)致了質(zhì)量流速和熱提取速率的降低，壓力變化幅度減小?？梢娭車紫督橘|(zhì)對壓裂區(qū)域傳導(dǎo)的熱量還是占有一定比例的。因此，在深部地?zé)衢_采模擬過程中應(yīng)該考慮上下基質(zhì)傳熱過程。

圖14 裂隙介質(zhì)裂隙度對溫度、壓力和熱提取效率的影響Fig.14 Effects of porosity of fractured media on temperature，pressure，and mass and heat extraction rate

通過參數(shù)的不確定分析可知，裂隙系統(tǒng)的滲透率、注入井/生產(chǎn)井壓力和注入溫度、井徑對深部地?zé)衢_采過程中的壓力和溫度影響較大，從而影響到熱的提取效率，應(yīng)重視這些因素對松遼盆地地?zé)衢_采過程的影響。而孔隙介質(zhì)的滲透率和孔隙度、裂隙介質(zhì)的孔隙度和巖石的熱傳導(dǎo)系數(shù)的影響并不明顯。

5 結(jié)論與展望

1）松遼盆地徐家圍子斷陷區(qū)地溫梯度較高，埋深4 100m地層溫度達到160℃，是良好的熱儲層；其巖性主要為砂巖，具有一定的孔隙度和滲透率（孔隙度在0.1左右，滲透率在10－17m2數(shù)量級上），但仍然需要采用水力壓裂來人工增加滲透性，以達到規(guī)?；_發(fā)的目的。

圖15 巖石熱傳導(dǎo)系數(shù)對溫度、壓力和熱提取效率的影響Fig.15 Effects of rock thermal conductivity on temperature，pressure，and mass and heat extraction rate

2）在采用定壓力開采地?zé)岬那闆r下，注入井和生產(chǎn)井區(qū)域溫度的差異使得水的流動性差別較大，生產(chǎn)井控制整個區(qū)域的壓力分布，壓力梯度在注入井區(qū)域較大，隨著開發(fā)的進行，注入井的影響逐漸增大；裂隙介質(zhì)中的流體溫度由注入井到生產(chǎn)井逐漸升高，并隨著開發(fā)的進行而降低，逐漸向生產(chǎn)井發(fā)展，其中的孔隙基質(zhì)溫度變化特征與流體相似，并且在溫度變化前緣存在明顯的梯度，即離裂隙越近，溫度降低越大。

3）不同埋深位置的模擬顯示，埋深大的溫度相對較高，水的流動性較強，質(zhì)量和熱提取速率較高，壓力和溫度變化均較大。

圖16 井徑對溫度、壓力和熱提取效率的影響Fig.16 Effects of well radius on temperature，pressure，and mass and heat extraction rate

4）裂隙系統(tǒng)的滲透率、注入井/生產(chǎn)井壓力和注入溫度、井徑對深部地?zé)衢_采過程中的壓力和溫度影響較大，從而影響到熱的提取速率，而孔隙基質(zhì)的滲透率和孔隙度、裂隙介質(zhì)的孔隙度和巖石的熱傳導(dǎo)系數(shù)的影響并不明顯。

本次研究結(jié)論是在特定參數(shù)條件下獲得的，并且只考慮了壓力場和溫度場的耦合效應(yīng)，沒有綜合考慮應(yīng)力場和化學(xué)場的影響，而實際工程條件下這四個場的互相影響是客觀存在、復(fù)雜和重要的，因此，還需在以后的研究中進一步深入。

圖17 注入井/生產(chǎn)井壓力對溫度、壓力和熱提取效率的影響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 壓裂區(qū)域上下基質(zhì)對溫度、壓力和熱提取效率的影響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］許天福，金光榮，岳高凡，等.地下多組分反應(yīng)溶質(zhì)運移數(shù)值模擬：地質(zhì)資源和環(huán)境研究的新方法［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，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，F(xiàn)rance［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］高瑞琪，蔡希源.松遼盆地油氣田形成條件與分布規(guī)律［M］.北京：石油工業(yè)出版社，1997.Gao Ruiqi，Cai Xiyuan.Field Formation Conditions and Distribution Rules in Songliao Basin［M］.Beijing：Petroleum Industry Press，1997.

［12］周慶華，馮子輝，門廣田.松遼盆地北部徐家圍子斷陷現(xiàn)今地溫特征及其與天然氣生成關(guān)系研究［J］.中國科學(xué)：D輯：地球科學(xué)，2007，37（增刊II）：177－188.Zhou Qinghua，F(xiàn)eng 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］任戰(zhàn)利.中國北方沉積盆地構(gòu)造熱演化史研究［M］.北京：石油工業(yè)出版社，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，F(xiàn)ujimitsu 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.