蒙特卡罗法在地热资源评价中的应用

韩 军， 刘 祜， 甄 珍， 朱鹏飞

(1. 中核集团铀资源勘查评价重点实验室，北京 100029； 2. 核工业北京地质研究院，北京 100029)

0 引 言

在热储概念模型基础上，根据地热条件和工作程度，参照现有地热标准(GB/T 11615—2010、DZ/T 40-85)，结合国内外地热资源估算最新方法得出地热资源量评价(计算)方法。地热资源评价方法主要有热储法(体积法)、端压法、解析法、比拟法、统计分析法、数值模拟法等，实际工作一般根据地热勘查工作阶段和掌握的地热参数来确定采用何种评价方法。

根据印尼B地热田的踏勘结果、地热工作程度和已有勘查评价资料，参照我国地热相关行业标准，采用基于热储法的蒙特卡罗(Monte Carlo)法初步评价了该地热田的可利用地热资源量，满足了开发方要求，并讨论了该方法的适用条件和结果，为下一步工作提供可靠依据。

1 基本地质条件和地热勘查程度

根据我国地热勘查规范《地热资源地质勘查规范》(GB/T 11615—2010)，研究区地热勘查工作网度和勘查手段均处于调查阶段，勘探程度较低。利用已有网度较稀疏的电磁法剖面和地热化探测量结果初步圈出热储范围与埋深，并根据少量水化学分析结果初步预测热储温度。由于工作程度低，采样点数量较少，且地热田及周边没有钻探控制，因此达不到准确预测程度。研究区地质和地热基本条件如下。

(1) 研究区为受构造控制的地热异常区，发育火山岩和火山机构，在出露温度较高的B、M两个温泉点(图1、图2)，出露泉点最高水温分别为83.0、44.3 ℃，均为变质岩和岩浆岩基岩裂隙型热储，初步预测浅层热储埋深为800～1 000 m。

图1 M-B温泉区域地质简图(据印度尼西亚能源和矿产资源部，2011)1-第四纪冲洪积物；2-第四纪粗碎屑沉积物；3-中新世花岗岩；4-千枚岩；5-板岩；6-古生代花岗岩；7-片岩；8-断层；9-地热点Fig. 1 Regional geological sketch map of M-B thermal spring area(after Indonesian Ministry of Energy and Mineral Resources, 2011)

(2) 两温泉点主要蚀变矿物为明矾石、蒙脱石、高岭土类黏土矿物，蚀变矿物属于低温—中高温(50～300 ℃)热液作用的产物(王联魁等，1977)，B温泉点蚀变矿物的分布范围和厚度较大。

(3) 根据水化学分析结果，B温泉点具深部热水来源，与流经地层发生的水岩交换程度较低；M温泉水与B温泉水的水化学类型、平衡程度差别较大，具不同热储来源，后者与浅层水混合程度较大。

(4) 地表出露热水温度最高达90～93 ℃，地热异常分布区总体面积>25 km2，预测浅部(1 km以浅)范围为10 km2，热储估算温度为134～220 ℃。

图2 M-B温泉区剖面示意图(据印度尼西亚能源和矿产资源部，2011)Qalp-第四纪冲洪积物；Qclsd-第四纪粗碎屑沉积物；Kf-中生代千枚岩；Ks-中生代板岩；Trs-中生代片岩；Tgs-中新世花岗闪长岩、闪长岩；Tgo-中新世花岗岩；Trg-古生代花岗岩、花岗片麻岩Fig. 2 Section sketch of M-B thermal spring area(after Indonesian Ministry of Energy and Mineral Resources, 2011)

(5) 根据地热异常推测深部热储中心位于B温泉附近，据物探电法测量切片反演结果，热储位于B、L区附近2组断裂延长线地段，即现有B、L区的北侧，由于该区物探测量网度稀疏，需要加密物探测量证实。

(6) 上述计算仅针对1 km以浅的浅层热储，热储面积为B、L区所在的小片区域，考虑到深部岩浆侵入带来的热源和提供的地热流体具有更高的温压值，因此估算的资源量较为保守。

(7) 根据大地构造环境、构造岩浆演化以及已有钻孔资料分析(印度尼西亚能源和矿产资源部，2011)，该区深部存在火山岩浆热源的裂隙型热储，B温泉周围深部存在>200 ℃的基岩裂隙型热储，具有有利的变质岩热储盖层以及由断裂构造沟通的地热流体循环、补给系统，预测在1 km以浅可揭遇180 ℃的热水，在>1 km深处具有深部第二热储，地热开发潜力显著。

2 基于热储法的蒙特卡罗法模拟

2.1 热储法基本原理

热储法普遍适用于各类、各阶段地热资源储量的评价(韩征等，2015；王彩会等，2015)，运用蒙特卡罗法时以该方法为基本公式进行模拟计算(Lee et al., 1980)。热储法将大部分浅层热储视为圈闭型，即热能补给充分，水补给少，基本为不可再生型地热资源。将地热资源量视为地热流体存储量，即热储及流体范围内的静态储量，不考虑侧向补给、越流补给，将热储层概化为平面上无限延伸、垂向上2层隔水层夹1个主要含水层的3层结构体系(王心义等，2002)。根据热储面积和热储有效厚度(评价基准面顶底面埋深、孔隙度、经孔隙度修正后的有效厚度)，确定热储几何形状、温度、孔隙度。将评价区分为若干子区，计算各子区热流体储量，分别计入权重后加和，得到评价区热流体总储量(汪集旸等，1993；闫晋龙等，2020)。计算公式为：

Q=Crρr(1-φ)V(T1-T0)+Cwρwqw(T1-T0)

(1)

qw=φV+μA(h-H)

(2)

式中：Q为地热储量，kJ；Cr、Cw分别为热储层岩石和地热流体的热导率，kJ/(kg·℃)；φ为热储岩石(层)孔隙度或裂隙率；ρr、ρw分别为热储层岩石和地热流体的密度，kg/m3；T1、T0分别为热储温度、尾水温度或恒温层温度，℃；qw为地热流体弹性热存储量，包括弹性储量和静态储量，m3；A为热储面积，m2；V为热储体积，m3；h为平均承压水头高度，m；H为平均热储顶板算起的水头高度(一般用开采后允许的最大降深代替，或根据地热井水文试验数据),m；μ为含水层弹性释水系数。

2.2 蒙特卡罗法统计分析应用

2.2.1 基本原理 蒙特卡罗法自1949年提出并在地热资源预测和评价中得到了应用与发展，采用热储法(体积法)进行统计分析(Metropolis et al.，1949；Kappelmeyer et al.,1974；White et al.，1975；Brook et al.,1978；Nathenson，1978)，即根据热储法计算主要参数的数值范围、不同参数在公式中的权重等，获得不同概率的评价结果，以最大值、最小值、平均值来表示热储法中的3个重要参数(温度、体积和热能量)(王心义等，2002)。

蒙特卡罗法由Nathenson(1978)和Sarmiento等(2007)发展成熟，其核心是为每一个系统的主要参量(热储温度、面积、厚度等)估算出3个数值，作为三角形概率密度的最小、最可能和最大的估计值。估值范围接近变量真实概率密度的估计值。可运用解析方法计算设定的三角形概率密度的平均值和标准偏差，从而得到目标地热田的温度、体积、能量平均值和标准差。该方法的优点是可给出待评价资源量的极值范围、取值置信度和概率分布，提供预测值范围和最乐观值、最保守值等数据(Rybach et al.，1986)。

在使用热储法进行地热资源评估时，一般将参数设为定值进行计算。不同勘查阶段、不同地热田获得的热储参数值通常为一个范围，即在一定范围内变换的随机变量。在进行资源量评估时，热储的厚度、面积、温度、孔隙度等主要参数均为范围值和可能值，其精度决定了资源量评估结果的精度。蒙特卡罗法主要基于热储法(体积法)设定随机变量和预测结果的性质、范围，以概率及统计学的理论方法为基础，根据待求随机问题的变化规律和物理现象自身的统计规律构造概率模型，然后通过统计实验使某些统计参量恰为待求问题的解(李同彪，2015)。

蒙特卡罗法在模拟计算地热资源量(预测发电量)时，一般通过设定随机过程反复生成时间序列，设定参数估计量和统计量(统计结果)，同时给出结果的概率分布模式(函数)。该方法可利用较多的热储资料，较好地解决计算结果的可信度问题(Nathenson，1978)。

模拟前应设计模型(模拟公式)，设定公式中的随机变量(变量)，提出问题(预测量)，设定变量和预测量的概率分布模型(分布函数)；然后给出各随机变量的值域和最可能取值，设计某种规则(模拟公式)，根据模型中各随机变量的分布类型，通过大数据量的模拟试验(随机抽样试验)对随机抽样数据进行数学计算、统计分析模拟试验，使预测量的解对应于模型中随机变量的某些特征(如概率、均值等)，得出的预测量通过概率分布进行显示(颜英军，2009；李同彪，2015)。

试验结果包括最小值、最大值、数学期望值和标准偏差等，生成概率分布曲线图以定义结果的类型(李同彪，2015)。

2.2.2 应用实例 目标地热田为隆起山地对流型地热资源，属于变质-火山岩基岩裂隙型热储。预测热储埋深在800～2 000 m之间，可视为圈闭型热储。基于热储法，将公式(1)变换以适于蒙特卡罗法的应用(颜英军，2009；韩征等，2015)。

Q储=Q容=XAHφ

(3)

式(3)中：Q储为地热流体存储量，m3；Q容为地热流体容积储量，m3；A为热储面积，m2；H为热储层厚度，m；φ为热储层孔隙度；X为可采系数，受热储层岩性、厚度、孔隙度等控制，浅层水的开采系数X取值一般在1%～2%之间，深层水X值稍大于浅层水(邱楠生等，2004)。

利用蒙特卡罗法模拟，设定95%置信度条件，给出100%、95%、75%、50%概率条件下的预测资源量(可发电量)。热储法模拟计算需要的参数(条件)有热储面积和厚度、热储参数、渗透率，以及基础岩石学、热力学、水力学参数等。

(1) 模拟公式及变量预测量定义。为研究地热田地热资源可利用量(发电量)，基于公式(3)设计蒙特卡罗法模拟计算的预设变量和预测变量，计算公式为：

(4)

式(4)中：Q为实际可发电量，MW；k为综合利用系数，根据地热采收率和发电效率综合而得；A为预测的热储面积，m2；H有为热储有效厚度，m；C综为体积法计算的综合热容比，kJ/(kg·℃)；Tpg为热储温度，℃；Tcut-off为根据地热发电利用确定的地热流体下限温度，℃，取120 ℃(印度尼西亚能源和矿产资源部，2011)；(30×7 000×3 600)为地热发电时间(30年发电期，每年利用时长7 000 h)，s;×1 000表示将发电量结果换算成MW。

表1 B地热田蒙特卡罗法模拟计算参数值

H有=αH

(5)

式(5)中:α为储厚比，参照基岩裂隙型热储的经验值，介于0.02～0.40之间；H为热储层所在地层总厚度，m。

C综=ρrCr(1-φ)+ρwCwφ

(6)

式(6)中：ρw为T温度时水的密度，kg/m3，估算热储温度在100～200 ℃之间，取定值920；ρr为T温度时岩石的密度，kg/m3，热储地层为变质岩和侵入岩，取定值2 500；φ为孔隙度，根据钻井实测值确定取值范围为0.04～0.19；Cw为水热容，kJ/(kg·℃)，取定值4.18；Cr为岩石热容，kJ/(kg·℃)，采用定值0.879。

模拟时定义了资源量计算主要参数的函数类型(三角函数、正太函数、对数函数等)，确定合理的取值(最大值、最小值、平均值、众数值等)范围，参考地热地质条件对模拟结果进行分析，保证结果具有较好的可靠性。

发电量Q根据公式(4)、(5)、(6)模拟计算而得，式中采用的变量、预测量及范围见表1。

(2) 模拟结果。使用Oracle Crystal Ball软件,设置300 001次随机模拟次数，计算B地热田蒙特卡罗法的预测结果，设置为95%的置信度区间，预计30年发电期，每年利用7 000 h，发电量单位为MW。具体结果如下。

① 100%概率条件下的预测值。图3显示了一种常用集中分布函数，属于正态分布类型(正态函数)，结果为对称分布。最可能发生值是图中显示的平均值59.14 MW，标准偏差为30.28。图3显示，接近平均值的数发生几率较高。

图3 100%概率条件下B地热田预测发电量Fig. 3 Predicted electric generating capacity of the B geothermal field under the probability of 100%

根据切比雪夫定理(王学仁，1982)，任何样本中，落在平均值的k个标准偏差范围内的观测值例(几率)至少为1-1/k2，其中k>1。当k=3时，几率>8/9，因此以平均值±3倍标准偏差为预测值的最大值，取值范围为0～149.98 MW。

② 99%概率条件下的预测值。图4显示了正态分布函数。最可能发生值是平均值59.14 MW，在该概率下预测的最小值为8.92 MW，最大值为157.56 MW。

图4 99%概率条件下B地热田预测发电量Fig. 4 Predicted electric generating capacity of the B geothermal field under the probability of 99%

③ 95%概率条件下的预测值。图5显示了一种正态分布函数。最可能发生值是平均值59.14 MW，该概率下预测的最小值为14.25 MW，最大值为129.55 MW。

图5 95%概率条件下B地热田预测发电量Fig. 5 Predicted electric generating capacity of the B geothermal field under the probability of 95%

④ 75%概率条件下的预测值。图6显示了正态分布函数。最可能发生值是平均值59.14 MW，在该概率下预测的最小值为26.18 MW，最大值为95.26 MW。

图6 75%概率条件下B地热田预测发电量Fig. 6 Predicted electric generating capacity of the B geothermal field under the probability of 75%

⑤ 50%概率条件下预测值。图7显示了正态分布函数。最可能发生值是平均值59.14 MW，该概率下预测的最小值为36.3 MW，最大值为77.2 MW。

图7 50%概率条件下B地热田预测发电量Fig. 7 Predicted electric generating capacity of the B geothermal field under the probability of 50%

⑥ 随机变量权重分析。B地热田预测模型敏感度图解(图8)显示，热储面积参数V2占比55.5%，热储温度参数V6占比37.6%，综合利用系数参数V8占比3.4%，有效热储厚度参数V3占比3.3%，综合热导系数参数V7占比0.2%。

敏感度图显示了模拟公式中随机变量参数对结果的影响程度。图8表明，热储面积、热储温度、利用系数、热储厚度、综合导热系数5项参数为模拟结果(发电量)的主要影响因子，且对结果的影响程度依次降低。

图8 B地热田预测模型敏感度图解Fig. 8 Sensitive chart of predicted model of the B geothermal field

3 讨 论

3.1 预测发电量最小值

图3—图7为蒙特卡罗法在不同概率下的结果和函数分布类型，根据正态函数性质，100%概率模式下的最小值>0，但未给出实际值。

根据地表热释放量法计算的结果为834.86 kW(印度尼西亚能源和矿产资源部，2011)，根据世界范围内的经验值(中国科学院地质研究所地热组，1978)，地面释放量一般比实际值至少差1～3个数量级(李同彪，2015)，以10倍对其修正后得出的值(8.348 6 MW)作为该方法可利用资源量的最小下限值。

运用体积法，设计发电期30年、年利用时长为7 000 h、SiO2的热储温度为161 ℃、回收温度为120 ℃、孔隙度为10%，估算资源量为27.00 MW。

两者对比，体积法更接近且小于实际预测的发电量，因此以27.00 MW为预测发电量下限值。

3.2 蒙特卡罗法预测值

蒙特卡罗法可用于地热开发各阶段资源量的预测(王心义等，2002；韩征等，2015)，对有多年动态监测资料的地热田的资源量估算更准确(颜英军，2009)，在估算资源量时可给出置信区间(李同彪，2015)。研究区尚无地热钻孔，有效的地热资源计算方法为体积法，但其计算参数均为定值，无法给出置信区间，因此采用已有资料和分析数据时均选用保守参数值来计算，可认为在目前地热调查阶段59.14 MW的地热发电量可信度高。随着勘查工作的深入，其资源量呈增加趋势。

4 结 论

(1) 此次地热资源量预测仅针对1 km以内的浅层热储，热储面积仅为B温泉所在的小片区域。

(2) 结合体积法固定参数计算和模拟试验结果分析，认为资源量计算的最可能值在置信度为95%时为59.14 MW，可能范围为27.00～149.98 MW，远景发电量为115.85 MW。

(3) 影响发电量预测值的主要参数随机率降低依次为预测热储面积、热储温度、利用系数、热储厚度、综合导热系数。

(4) 根据大地构造环境、构造岩浆演化及邻区钻孔资料分析，研究区深部存在火山岩浆热源的裂隙型热储，特别是B温泉区深部存在>200 ℃的基岩裂隙型热储，具有有利的变质岩热储盖层和由断裂构造沟通的地热流体循环、补给系统，在深度>1 km的深部有第二热储，具有较大的地热资源开发潜力。