蒙特卡羅法在地熱資源評價中的應用

韓 軍， 劉 祜， 甄 珍， 朱鵬飛

(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的深部有第二熱儲，具有較大的地熱資源開發潛力。