2. 吉林大学地下水资源与环境教育部重点实验室, 长春 130021;
3. 大庆油田井下作业分公司, 黑龙江大庆 163453
2. Key Lab of Groundwater Resource and Environment, Ministry of Education, Jilin University, Changchun 130021, China;
3. Downhole Operation Company of Daqing Oilfield, Daqing 163453, Heilongjiang, China
0 前言
大庆油田的常年开采不仅造成石油产量逐渐减少，而且开采所用的化石燃料燃烧也产生了大量温室气体;因此，需要寻找一种清洁能源来代替化石燃料。油井测温数据显示，大庆徐家围子地区地温梯度为3.8~6.5 ℃/hm，2 km以下地温梯度稳定在3.8~4.0 ℃/hm，4 km深处温度约160 ℃[1]，具有非常高的深部地热能开采价值。目前开采深部地热(干热岩)的技术称为增强型地热系统(enhanced geothermal system,EGS)[2]，即通过对天然低渗地层进行人工改造，使之形成具有一定规模和足够渗透能力的天然热储层，注入冷水通过热储层加热，抽出热水用来发电和供热。
储层改造是EGS工程中非常重要的一个部分，主要包括化学刺激、热刺激和水力压裂，以水力压裂为主。水力压裂根据破裂机理，可以分为张开型压裂和剪切型压裂；根据压裂液类型，可以分为支撑剂型和清水型[3]。张开型压裂的应用非常广泛，几乎对地层条件没有要求；剪切型压裂应在天然裂缝发育的储层，配合热刺激和化学刺激进行。目前国际上EGS工程的开发主要是选择在天然裂隙发育的地层、采用清水剪切型压裂为主的储层改造方式，而传统的支撑剂型压裂由于改造体积有限和换热效果差的原因不受重视。
然而，从全球EGS发展历程来看，EGS的开发并没有完全固定的模式，而是具有“因地制宜”性，因此越来越多的研究倾向于根据地方需求进行EGS开采。Van der Hoorn等[4]基于水力压裂模拟研究了荷兰2个石灰岩EGS工程的换热情况。Hannes Hofmann等[3]针对加拿大Alberta油田深部低温(100 ℃)花岗岩储层是否可以为石油工程提供50~60 ℃热水的问题，进行了水力压裂和换热研究。Zeng等[5]采用水平井技术对Desert Peak场地进行了水热产出预测并进行了参数敏感性分析。可见，只要EGS工程通过储层改造能够达到工程经济性应用的目的就很有价值。对于天然裂隙发育的储层，应优先考虑采用清水剪切型压裂(如Desert Peak)形成裂隙网络产出高温热水来发电；对致密储层，可以考虑采用张开型压裂(像芬顿山或石油压裂)形成裂缝型储层产出中低温热水来供热。因此，本文基于大庆莺深2井的测井资料，对其深部存在的各种储层条件进行了划分，分别提出了相应的储层改造方案并建立了三维压裂模型，基于不同储层改造结果，利用水热耦合软件TOUGH2-EOS1进行了换热分析，得出了不同储层的干热岩开采潜力。
1 地质背景 1.1 地质构造莺深2井位于双城断陷区的莺深井区，该区受4条近NNE走向的“S”形断裂控制，呈现“两凹一凸”的构造格局[6](图 1)，表现为双重结构的地质特征，分为上部凹陷期和下部断陷期。断陷期自下而上主要包括火石岭组、沙河子组、营城组[6]。该区营城组发育较好、分布广，主要岩性为流纹岩、凝灰岩、安山岩、砂泥岩。沙河子组下段下部为泥岩夹砾岩、粉砂质泥岩，中上部为泥岩；上段下部为厚层泥岩夹砾岩，上部为厚层泥岩与砾岩互层。火石岭组上部和下部是中基性火山岩，中部是碎屑岩夹煤层。
 |
| 图 1 双城断陷勘探构造图Fig.1 Structure map of Shuangcheng fault depression |
根据从油田获得的测温数据及参考周庆华等[7]研究成果，绘制了徐家围子地区不同深度的地温分布图(图 2)，得出现今徐家围子地区地温梯度约为4 ℃/hm；3 km深处地温108~136 ℃，平均121 ℃；4 km深处144~180 ℃，平均162 ℃。
 |
| 图 2 徐家围子不同深度的地温分布图Fig.2 Temperature distribution map at different depths in Xujiaweiz area |
吉林大学物探学院对莺深2井(YS2)进行了大地电磁测量，解译结果(图 3)显示：2 km浅层A区低阻层对应莺山凹陷内的晚期沉积地层；4 km深处B区左侧的高阻区等值线形态较好地对应了凹陷区东侧的青山凸起；深部出现了局部的低阻异常区C区，推断可能为深部地热来源。
 |
| 图 3 莺深2井MT解译图Fig.3 MT interpretation results of YS 2 well |
干热岩要求储层温度大于150 ℃，对应该地区4 km及以下的地层，该层位为营城组巨厚火山岩(3 784.0~5 020.0 m)[7]；因此，从深度、温度和热源三方面来说，营城组火山岩储层是该区目前最适合进行干热岩开采的储层。
1.3 营城组储层特征通过对莺深2井营城组地层进行成像测井，得到了3 865.0~4 493.3 m深度内的裂缝发育情况(表 1)。Soultz场地干热岩压裂储层的裂隙节点密度为3.0个/m，因此我们以此为标准来选择该井适合进行干热岩压裂开发的裂隙发育段。由表 1可看出，裂隙发育集中在3 883.7~4 064.8 m和4 214.3~4 323.4 m两段内(平均裂隙节点密度大于3.0个/m)，其余地层岩体致密。对于致密岩体段，由于裂隙极其不发育，只能采用张开型压裂方式对其进行改造；对天然裂隙发育的地层段，应优先采用剪切型压裂，配合化学刺激和热刺激对其进行改造。
| 井段/ m | 裂缝长度/ m | 裂缝水动力 宽度/mm | 节点密度/ (个/m) |
| 3 865.0~3 867.4 | 1.5 | 0.9 | 2.1 |
| 3 883.7~3 909.2 | 3.4 | 0.5 | 4.3 |
| 3 914.2~3 954.6 | 3.4 | 0.3 | 4.4 |
| 3 957.0~3 971.7 | 1.8 | 0.3 | 2.4 |
| 3 973.9~4 023.2 | 2.4 | 0.4 | 3.1 |
| 4 029.0~4 030.0 | 1.4 | 0.5 | 2.0 |
| 4 037.3~4 045.1 | 2.4 | 0.5 | 3.0 |
| 4 059.1~4 064.8 | 2.6 | 0.4 | 3.4 |
| 4 076.4~4 081.8 | 1.3 | 0.3 | 1.8 |
| 4 087.8~4 130.0 | 2.3 | 0.2 | 3.1 |
| 4 144.6~4 150.0 | 1.9 | 0.6 | 2.7 |
| 4 178.3~4 206.2 | 1.7 | 0.2 | 2.2 |
| 4 214.3~4 234.1 | 2.4 | 0.7 | 3.1 |
| 4 245.5~4 258.2 | 4.2 | 0.2 | 5.5 |
| 4 269.9~4 283.0 | 2.3 | 0.3 | 3.1 |
| 4 287.7~4 307.9 | 3.1 | 0.2 | 4.1 |
| 4 316.2~4 323.4 | 2.4 | 0.3 | 3.3 |
| 4 335.7~4 346.0 | 1.6 | 0.2 | 2.1 |
| 4 351.6~4 384.5 | 1.8 | 0.2 | 2.5 |
| 4 390.8~4 401.0 | 1.5 | 0.1 | 2.2 |
| 4 410.6~4 415.6 | 1.5 | 0.2 | 2.0 |
| 4 443.0~4 446.8 | 1.9 | 0.4 | 2.4 |
| 4 459.9~4 493.3 | 3.0 | 0.2 | 4.0 |
本次张开型压裂采用认可度非常高的商业化石油压裂软件STIMPLAN[8]来进行压裂建模及分析。张延军等[9]利用该软件模拟了EGS的水力压裂，证明了该软件在EGS压裂中具有适用性。
2.1 致密储层压裂模拟 2.1.1 建立压裂模型根据地应力变化情况可将目标储层分为两部分：遮挡层(3 870.0~3 925.8 m)和无遮挡层(3 925.8~3 950.0 m)(图 4)。对于遮挡层，最小主应力为53.0 MPa；上部3 870.2 m处为82.7 MPa，下部3 925.8 m处为65.0 MPa；上下应力遮挡较好，但应力波动较大。从密度曲线也可看出遮挡层密度较邻层低，物性较差，压裂过程中会很好地抑制裂缝高度(纵向)的发育，而使裂缝长度更长。对于无遮挡层，地应力几乎随深度成线性增长(0.025 MPa/m)，压裂过程中裂缝会向上延伸，高度变大，长度变小。由于这两种储层条件对裂缝几何形状影响很大，因此我们分别进行了压裂和热产出分析。相关储层参数见表 2。
 |
| 注：①in(英寸)为非法定计量单位，1 in=0.025 4 m,下同。图 4 目标储层地应力曲线Fig.4 In-stress curve in target formation |
| 井段/ m | 杨氏模量/ GPa | 泊松比 | 孔隙度/ % | 渗透率/ mD① |
| 3 871.6~3 873.6 | 43.2 | 0.314 | 5.8 | 0.10 |
| 3 879.4~3 910.6 | 54.8 | 0.188 | 8.3 | 0.32 |
| 3 919.8~3 953.4 | 59.5 | 0.192 | 2.9 | 0.01 |
| 3 968.4~4 015.4 | 57.3 | 0.175 | 2.6 | 0.01 |
| 4 089.8~4 132.6 | 62.4 | 0.220 | 1.3 | 0.01 |
| 4 138.2~4 140.2 | 55.8 | 0.201 | 4.1 | 0.06 |
| 4 213.6~4 324.0 | 60.1 | 0.126 | 1.8 | 0.01 |
| 注:①mD(毫达西)为非法定计量单位，1 mD=0.987×10-15 m2,下同。 | ||||
根据上述储层参数和大庆油田压裂经验，我们建立了储层水力压裂模型，模型参数见表 3。采用阶梯式逐级提高支撑剂浓度的方式进行压裂，每注入200 m3压裂液提高一次支撑剂质量浓度。
| 参数 | 值 |
| 层位 | 遮挡层和无遮挡层 |
| 支撑剂 | HSP 16/20 |
| 压裂液 | 50#X-Link |
| 排量/(m3/min) | 6.5 |
| 压裂液体积/m3 | 2 000 |
| 前置液比/% | 30 |
| 支撑剂质量浓度/ (kg/m3) | 87、121、173、242、 277、311、382、417 |
| 模拟时间/min | 307.3 |
图 5是遮挡层压裂裂缝有效宽度和导流能力图。由于良好的应力遮挡，裂缝半长几乎是高度的6倍，裂缝长约370 m，高约65 m，开度约2.5 mm，导流能力约1 440 mD·m。
 |
| 图 5 遮挡层形成的裂缝有效宽度(a)和导流能力(b)图Fig.5 Fracture aperture (a) and conductivity(b) map in confined formation |
图 6是无遮挡层压裂裂缝有效宽度和导流能力图。由于地应力随深度呈梯度增加，上下遮挡较弱，裂缝上下均可延伸，裂缝半长相当于高度的1/2。裂缝长约117 m，高约240 m，开度约2.1 mm，导流能力约1 120 mD·m。
 |
| 图 6 无遮挡层形成的裂缝有效宽度(a)和导流能力(b)图Fig.6 Fracture aperture (a) and conductivity (b) map in unconfined formation |
根据上述压裂结果，将生成的压裂储层概化成遮挡层和无遮挡层2种情况(表 4)。由于压裂模型各向同性，因此产生的裂缝关于井筒对称，概化的储层长度相当于2个裂缝半长。
| 地层 | 深度/ m | 长度/ m | 高度/ m | 宽度/ mm | 渗透率/ D | 最小水平地 应力/MPa |
| 遮挡层 | 3 872~3 922 | 600 | 50 | 2.5 | 576 | 53 |
| 无遮挡层 | 4 070~4 270 | 200 | 200 | 2.0 | 533 | 62 |
采用三水平井的开采模式，中间注两边抽且仅注水井进行压裂改造，水平井的水平段垂直裂缝面穿过。根据单缝换热情况来决定压裂次数和水平井水平段的长度，开采模式见图 7。
 |
| 图 7 三水平井模式示意图Fig.7 Schematic plot of 3-horizontal well systems |
本次采用美国劳伦斯伯克利国家重点实验室开发的多相流多组分软件TOUGH2的EOS1模拟器[17]，该模拟器对地下水热耦合模拟的准确性已经被很多人证实。根据表 4的参数分别建立遮挡层和无遮挡层的水热耦合模型。由于石油驱油加工和供暖均需要50 ℃以上的热水，因此要求EGS工程运行20 a产出温度不低于50 ℃。采用大庆地区地表水作为注入水，大庆年平均气温约5 ℃，假定注入水达到储层深度时被加热达到10 ℃，为简化起见，不考虑滤失。
2.2.2 换热模拟结果及讨论1)生产流速分析
图 8显示了遮挡层和无遮挡层注入压力随注入流速增加的变化规律。可以看出，注入压力随注入速率几乎呈线性增长。因最小水平地应力限制(表 4)，遮挡层和无遮挡层最大允许注入速率分别为6.7 kg/s和11.6 kg/s。无遮挡层注入压力上升速度比遮挡层的慢很多;原因是其注入井和生产井间距只有100 m，注入水很快被抽水井抽出，流动阻抗小。
 |
| 图 8 两种地层条件注入压力随注入速率的变化Fig.8 Injection pressure variation with the increase of injection flow rate in two formations |
图 9显示了20 a间遮挡层和无遮挡层出水温度随注入流速增加的变化规律。由于温度限制，可看出当注入速率达到4 kg/s时，两种储层第20年的出水温度均降到50 ℃，且遮挡层温降要比无遮挡层稍慢一些，随运行时间越长这种差距越来越小。
 |
| 图 9 两种地层条件出水温度随注入速率的变化Fig.9 Production water temeprature variation with the increase of injection flow rate in two formations |
2)产能分析
图 10是20 a间遮挡层和无遮挡层注入速率为1~4 kg/s的产能情况。可以看出，尽管随注入速率增加产出水温度降低(图 9)，但产热量反而增大，且遮挡层的产热量要高于无遮挡层。20 a间、4 kg/s的注入速率下，遮挡层和无遮挡层产热量分别是700~1 070 kJ/s和700~970 kJ/s。
 |
| 图 10 20 a间两种地层产出热量随注入速率的变化Fig.10 Production thermal energy variation with different injection flow rates in two formations |
3)水平井裂缝条数
遮挡层和无遮挡层的单缝换热在y方向的热传导影响范围分别是350 m(图 11a)和251 m(图 11b)。无遮挡层井间距较小，对围岩的热影响范围也较小。若采用水平井水平段2 000 m，两种储层可分别压裂5条和8条相互无热干扰的裂缝，相应的注入速率变为20 kg/s和32 kg/s，产热量变为3 500~5 350 kJ/s和5 600~7 760 kJ/s。可见，采用水平井技术时，选无遮挡层进行干热岩开采要好于遮挡层。
 |
| 图 11 遮挡层(a)和无遮挡层(b)在y方向热传导影响范围Fig.11 Heat conduction range in y direction of confined (a) and unconfined (b) formation |
从表 1可知，裂隙发育集中在3 883.7~4 064.8 m(I段)和4 214.3~4 323.4 m(II段)。虽然这两层裂隙较其他井段很发育，但其孔隙度(1.2%~8.3%)、渗透率(0.01~0.32 mD)依旧很低;因此,想要依靠该层的天然渗透能力循环水流提取地热能是不可能或不经济的，必须进行储层改造。
目前国际上对天然裂隙储层采用分阶段清水剪切压裂的方法进行改造。第一阶段先以低速率向地层泵入清水，随着注入速率的提高观察地表压力的变化。如果地表压力一直上升，则在第二阶段可以选择化学刺激注入酸液使裂隙扩大并连通，增大注入率。通过上述刺激如果储层流动阻抗仍大于0.1 MPa/(kg/s)[10]，则第三阶段采用张开型清水压裂进行改造，然后再重复进行第一、二阶段直到达到预定要求，典型的EGS工程代表比如Soultz和Desert Peak。因此，莺深2井的天然裂隙发育储层(I段和II段)也应采用这种方法来进行改造。
由于对远井的天然裂隙的密度、长度和展布特征并不了解，而剪切型压裂形成的储层跟天然裂隙网络的分布情况密切相关，因此无法采用数值模拟的手段来模拟剪切压裂后储层情况。然而，因压裂方法均是采用反复试探的手段使储层最终流动阻抗小于0.1 MPa/(kg/s)，因此我们将已有的EGS工程(Soultz和Desert Peak)与莺深2井进行对比，推断出压裂可能形成的储层情况(表 5)。结果显示，莺深2井经过相同改造后，在流动阻抗达到0.1 MPa/(kg/s)、注入速率为30 kg/s时，地表压力应不超过3 MPa，得到的EGS储层水平长2 000 m，水平宽200 m，垂直高181 m/109 m(I段/II段)，刺激体积达7.0×107 m3/5.6×107 m3(I段/II段)。采用的压裂及开采模式是中间注两边抽的三垂直井模式，井间距为500 m。
| Soultz | Desert Peak | 莺深2井 | |
| 深度/m | 5 000.0 | 1 300.0 | 4 000.0 |
| 岩性 | 花岗岩 | 流纹岩 | 流纹岩 |
| 天然裂隙层厚/m | - | 103 | 181/109 (I段/II段) |
| 天然裂隙节点密度/(个/m) | 3.2 | - | 3.3/3.9 (I段/II段) |
| 天然渗透率/mD | 0.01~0.10 | 0.10 | 0.30/0.01 |
| 刺激后流动阻抗/(MPa/(kg/s)) | 0.2 | 0.1 | 0.1 |
| 注入速率/(kg/s) | 30.0 | 18.0~54.3 | 30.0* |
| 地表压力/MPa | 7~8 | 3 | 3* |
| 刺激后水平长/m | 2 000 | 2 500 | 2 000* |
| 刺激后水平宽/m | 750 | 200或400 | 200* |
| 刺激后垂直高/m | 2 000 | - | 181/109(I段/II段) |
| 改造体积/m3 | 3×109 | - | 7.0×107/5.6×107(I段/II段) |
| 井间距/m | 500 | 1 000 | 500 |
| 注：“*”是与Soultz和Desert Peak对比得到的推断值(取二者的保守值)；Soultz和Desert Peak的参数来自文献[5,10,11]。 | |||
根据表 5，我们建立了剪切型压裂换热模型。注水段和抽水段与整个EGS储层厚度相同，I段和II段同时开采。本次模拟将EGS热储层等效为多孔介质储层，Pruess等[12]和Zeng等[5]等认为这种处理下得出的产能是可靠的，模拟参数见表 6。
| 参数 | 值 |
| EGS密度/(kg/m3) | 2 600 |
| EGS孔隙度 | 0.02 |
| EGS导热系数/(W/(m·K)) | 3.1 |
| EGS比热容/J/(kg·K) | 1 100 |
| 初始孔压/MPa | 采用公式P=43.23-10 000z计算 |
| 初始地层温度/℃ | 160 |
| 注水井速率/(kg/s) | 30 |
| 注水井地表压力/MPa | <3 |
| 注水井温度/℃ | 10 |
| 注：z为深度，单位为m。 | |
由3.1节可知，改造后储层渗透率应保证在20 a运行期间，注入速率为30 kg/s时注入压力升高不超过3 MPa。以注水井底部单元为观察点，初始平均孔压43 MPa，在20 a运行中，则不应超过46 MPa。通过不断调整EGS热储层渗透率，发现当渗透率为0.1 D时，注入压力在20 a期间刚好上升了3 MPa(图 12)。
随后我们进行了20 a出水温度和产能模拟，结果见图 13。可以看出:该储层出水温度比致密储层要高很多，在工程运行前10 a间，出水温度和产出热量都非常稳定，从第12年开始逐渐下降;20 a出水温度降低了10 ℃，产出热量为17 000~18 500 kJ/s。
 |
| 图 12 第20年储层孔压变化图Fig.12 Pore pressure value in 20th year |
 |
| 图 13 20 a水温度和产出热量变化Fig.13 Production water temperature and heat varation during 20 years |
1)营城组储层厚达1 236.0 m(3 784.0~5 020.0 m)，对应的温度很高(150~200 ℃)，是该区目前最适合进行干热岩开采的储层。
2)营城组目标储层根据天然裂隙发育程度，可以分为致密地层和天然裂隙发育储层；根据地应力遮挡情况，可以分为遮挡层和无遮挡层。
3)遮挡层和无遮挡层采用传统支撑剂型压裂改造后，在采用三水平井(一注两抽)(水平段2 000 m)的开采模式下，最大注入速率分别为20 kg/s和32 kg/s，20 a的产能分别为3 500~5 350 kJ/s和5 600~7 760 kJ/s。
4)天然裂隙储层经过清水剪切型改造后，20 a在注入速率30 kg/s下，产能达17 000~18 500 kJ/s。
5)虽然致密层产出热量远低于天然裂隙，但可以看出采用传统压裂和水平井技术后，使得致密储层也具有了开采干热岩的潜力。
| [1] | 雷宏武, 金光荣,李佳琦,等.松辽盆地增强型地热系统地热能开发热-水动力耦合过程[J]. 吉林大学学报(地球科学版), 2014,44(5):1633-1646. Lei Hongwu, Jin Guangrong, Lijiaqi, et al. Coupled Thermal-Hydrodynamic Processes for Geothermal Energy Exploitation in Enhanced Geothermal System at Songliao Basin, China[J]. Journal of Jilin University (Earth Science Edition), 2014,44(5):1633-1646. |
| [2] | 许天福,张延军,曾昭发,等. 增强型地热系统(干热岩)开发技术进展[J].科技导报,2012,30(32):42-45. Xu Tianfu, Zhang Yanjun, Zeng Zhaofa, et al. Technology Progress in an Enhanced Geothermal System (Hot Dry Rock)[J]. Science & Technology Review, 2012,30(32):42-45. |
| [3] | Hofmann H, Babadagli T, Zimmermann G. Hot Water Generation for Oil Sands Processing from Enhanced Geothermal Systems:Process Simulation for Different Hydraulic Fracturing Scenarios[J]. Applied Energy,2014,113:524-547. |
| [4] | Van der Hoorn K, Heijnen L J, Gankema M E, et al. Hydraulic Fracturing in Limestone:A Case Study of Two EGS Projects in The Netherlands[C]//Proceedings of the 37th Workshop on Geothermal Reservoir Engineering.Stanford:Stanford University,2012. |
| [5] | Zeng Y C, Wu N Y, Su Z, et al. Numerical Simulation of Heat Production Potential from Hot Dry Rock by Water Circulating Through a Novel Single Vertical Fracture at Desert Peak Geothermal Field[J]. Energy,2013,63:268-282. |
| [6] | Huang X X, Zhu J L, Niu C K, et al. Heat Extraction and Power Production Forecast of A Prospective Enhanced Geothermal System Site in Songliao Basin, China[J]. Energy,2014,75:360-370. |
| [7] | 周庆华,冯子辉,门广田. 松辽盆地北部徐家围子断陷现今地温特征及其与天然气生成关系研究[J].中国科学:D辑:地球科学,2007,37(增刊2):177-188. Zhou Qinghua, Feng Zihui, Men Guangtian. Research on Relationship Between Natural Gas Generation and Geothermal Characteristics in Xujiaweizi Depression Area in Northern Songliao Basin[J]. Science China:Series D:Earth Sciences, 2007,37(Sup.2):177-188. |
| [8] | Nolte K G, Smith M. STIMPLAN User's Guide[R]. Tulsa:NSI Technologies, Inc, 2011. |
| [9] | Zhang Y J,Li Z W,Guo L L,et al.Electricity Generation from Enhanced Geothermal Systems by Oilfield Produced Water Circulating Through Reservoir Stimulated by Staged Fracturing Technology for Horizontal Wells:A Case Study in Xujiaweizi Area in Daqing Oilfield, China[J]. Energy, 2014,78:788-805. |
| [10] | Chabora E, Zemach E, Spielman P, et al. Hydraulic Stimulation of Well 27-15, Desert Peak Geothermal Field, Nevada, USA[C]//Proceedings, Thirty-Seventh Workshop on Geothermal Reservoir Engineering. Stanford:Stanford University, 2012:SGP-TR-194. |
| [11] | Genter A, Fritsch D, Cuenot N, et al. Overview of the Current Activities of the European EGS Soultz Project:From Exploration to Electricity Production[C]// Proceedings, Thirty-Fourth Workshop on Geothermal Reservoir Engineering.Stanford:Stanford University, 2009:SGP-TR-187. |
| [12] | Pruess K, Oldenburg C, Moridis G. TOUGH2 User's Guide Version 2.0[R]. Berkeley:Lawrence Berkeley National Laboratory, 1999. |