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

(1. 中核集團(tuán)鈾資源勘查評價重點實驗室，北京 100029； 2. 核工業(yè)北京地質(zhì)研究院，北京 100029)

0 引 言

在熱儲概念模型基礎(chǔ)上，根據(jù)地?zé)釛l件和工作程度，參照現(xiàn)有地?zé)針?biāo)準(zhǔn)(GB/T 11615—2010、DZ/T 40-85)，結(jié)合國內(nèi)外地?zé)豳Y源估算最新方法得出地?zé)豳Y源量評價(計算)方法。地?zé)豳Y源評價方法主要有熱儲法(體積法)、端壓法、解析法、比擬法、統(tǒng)計分析法、數(shù)值模擬法等，實際工作一般根據(jù)地?zé)峥辈楣ぷ麟A段和掌握的地?zé)釁?shù)來確定采用何種評價方法。

根據(jù)印尼B地?zé)崽锏奶た苯Y(jié)果、地?zé)峁ぷ鞒潭群鸵延锌辈樵u價資料，參照我國地?zé)嵯嚓P(guān)行業(yè)標(biāo)準(zhǔn)，采用基于熱儲法的蒙特卡羅(Monte Carlo)法初步評價了該地?zé)崽锏目衫玫責(zé)豳Y源量，滿足了開發(fā)方要求，并討論了該方法的適用條件和結(jié)果，為下一步工作提供可靠依據(jù)。

1 基本地質(zhì)條件和地?zé)峥辈槌潭?/h2>

根據(jù)我國地?zé)峥辈橐?guī)范《地?zé)豳Y源地質(zhì)勘查規(guī)范》(GB/T 11615—2010)，研究區(qū)地?zé)峥辈楣ぷ骶W(wǎng)度和勘查手段均處于調(diào)查階段，勘探程度較低。利用已有網(wǎng)度較稀疏的電磁法剖面和地?zé)峄綔y量結(jié)果初步圈出熱儲范圍與埋深，并根據(jù)少量水化學(xué)分析結(jié)果初步預(yù)測熱儲溫度。由于工作程度低，采樣點數(shù)量較少，且地?zé)崽锛爸苓厸]有鉆探控制，因此達(dá)不到準(zhǔn)確預(yù)測程度。研究區(qū)地質(zhì)和地?zé)峄緱l件如下。

(1) 研究區(qū)為受構(gòu)造控制的地?zé)岙惓^(qū)，發(fā)育火山巖和火山機構(gòu)，在出露溫度較高的B、M兩個溫泉點(圖1、圖2)，出露泉點最高水溫分別為83.0、44.3 ℃，均為變質(zhì)巖和巖漿巖基巖裂隙型熱儲，初步預(yù)測淺層熱儲埋深為800～1 000 m。

圖1 M-B溫泉區(qū)域地質(zhì)簡圖(據(jù)印度尼西亞能源和礦產(chǎn)資源部，2011)1-第四紀(jì)沖洪積物；2-第四紀(jì)粗碎屑沉積物；3-中新世花崗巖；4-千枚巖；5-板巖；6-古生代花崗巖；7-片巖；8-斷層；9-地?zé)狳cFig. 1 Regional geological sketch map of M-B thermal spring area(after Indonesian Ministry of Energy and Mineral Resources, 2011)

(2) 兩溫泉點主要蝕變礦物為明礬石、蒙脫石、高嶺土類黏土礦物，蝕變礦物屬于低溫—中高溫(50～300 ℃)熱液作用的產(chǎn)物(王聯(lián)魁等，1977)，B溫泉點蝕變礦物的分布范圍和厚度較大。

(3) 根據(jù)水化學(xué)分析結(jié)果，B溫泉點具深部熱水來源，與流經(jīng)地層發(fā)生的水巖交換程度較低；M溫泉水與B溫泉水的水化學(xué)類型、平衡程度差別較大，具不同熱儲來源，后者與淺層水混合程度較大。

(4) 地表出露熱水溫度最高達(dá)90～93 ℃，地?zé)岙惓７植紖^(qū)總體面積>25 km2，預(yù)測淺部(1 km以淺)范圍為10 km2，熱儲估算溫度為134～220 ℃。

圖2 M-B溫泉區(qū)剖面示意圖(據(jù)印度尼西亞能源和礦產(chǎn)資源部，2011)Qalp-第四紀(jì)沖洪積物；Qclsd-第四紀(jì)粗碎屑沉積物；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) 根據(jù)地?zé)岙惓Ｍ茰y深部熱儲中心位于B溫泉附近，據(jù)物探電法測量切片反演結(jié)果，熱儲位于B、L區(qū)附近2組斷裂延長線地段，即現(xiàn)有B、L區(qū)的北側(cè)，由于該區(qū)物探測量網(wǎng)度稀疏，需要加密物探測量證實。

(6) 上述計算僅針對1 km以淺的淺層熱儲，熱儲面積為B、L區(qū)所在的小片區(qū)域，考慮到深部巖漿侵入帶來的熱源和提供的地?zé)崃黧w具有更高的溫壓值，因此估算的資源量較為保守。

(7) 根據(jù)大地構(gòu)造環(huán)境、構(gòu)造巖漿演化以及已有鉆孔資料分析(印度尼西亞能源和礦產(chǎn)資源部，2011)，該區(qū)深部存在火山巖漿熱源的裂隙型熱儲，B溫泉周圍深部存在>200 ℃的基巖裂隙型熱儲，具有有利的變質(zhì)巖熱儲蓋層以及由斷裂構(gòu)造溝通的地?zé)崃黧w循環(huán)、補給系統(tǒng)，預(yù)測在1 km以淺可揭遇180 ℃的熱水，在>1 km深處具有深部第二熱儲，地?zé)衢_發(fā)潛力顯著。

2 基于熱儲法的蒙特卡羅法模擬

2.1 熱儲法基本原理

熱儲法普遍適用于各類、各階段地?zé)豳Y源儲量的評價(韓征等，2015；王彩會等，2015)，運用蒙特卡羅法時以該方法為基本公式進(jìn)行模擬計算(Lee et al., 1980)。熱儲法將大部分淺層熱儲視為圈閉型，即熱能補給充分，水補給少，基本為不可再生型地?zé)豳Y源。將地?zé)豳Y源量視為地?zé)崃黧w存儲量，即熱儲及流體范圍內(nèi)的靜態(tài)儲量，不考慮側(cè)向補給、越流補給，將熱儲層概化為平面上無限延伸、垂向上2層隔水層夾1個主要含水層的3層結(jié)構(gòu)體系(王心義等，2002)。根據(jù)熱儲面積和熱儲有效厚度(評價基準(zhǔn)面頂?shù)酌媛裆睢⒖紫抖?、?jīng)孔隙度修正后的有效厚度)，確定熱儲幾何形狀、溫度、孔隙度。將評價區(qū)分為若干子區(qū)，計算各子區(qū)熱流體儲量，分別計入權(quán)重后加和，得到評價區(qū)熱流體總儲量(汪集旸等，1993；閆晉龍等，2020)。計算公式為：

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

(1)

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

(2)

式中：Q為地?zé)醿α?，kJ；Cr、Cw分別為熱儲層巖石和地?zé)崃黧w的熱導(dǎo)率，kJ/(kg·℃)；φ為熱儲巖石(層)孔隙度或裂隙率；ρr、ρw分別為熱儲層巖石和地?zé)崃黧w的密度，kg/m3；T1、T0分別為熱儲溫度、尾水溫度或恒溫層溫度，℃；qw為地?zé)崃黧w彈性熱存儲量，包括彈性儲量和靜態(tài)儲量，m3；A為熱儲面積，m2；V為熱儲體積，m3；h為平均承壓水頭高度，m；H為平均熱儲頂板算起的水頭高度(一般用開采后允許的最大降深代替，或根據(jù)地?zé)峋脑囼灁?shù)據(jù)),m；μ為含水層彈性釋水系數(shù)。

2.2 蒙特卡羅法統(tǒng)計分析應(yīng)用

2.2.1 基本原理 蒙特卡羅法自1949年提出并在地?zé)豳Y源預(yù)測和評價中得到了應(yīng)用與發(fā)展，采用熱儲法(體積法)進(jìn)行統(tǒng)計分析(Metropolis et al.，1949；Kappelmeyer et al.,1974；White et al.，1975；Brook et al.,1978；Nathenson，1978)，即根據(jù)熱儲法計算主要參數(shù)的數(shù)值范圍、不同參數(shù)在公式中的權(quán)重等，獲得不同概率的評價結(jié)果，以最大值、最小值、平均值來表示熱儲法中的3個重要參數(shù)(溫度、體積和熱能量)(王心義等，2002)。

蒙特卡羅法由Nathenson(1978)和Sarmiento等(2007)發(fā)展成熟，其核心是為每一個系統(tǒng)的主要參量(熱儲溫度、面積、厚度等)估算出3個數(shù)值，作為三角形概率密度的最小、最可能和最大的估計值。估值范圍接近變量真實概率密度的估計值?？蛇\用解析方法計算設(shè)定的三角形概率密度的平均值和標(biāo)準(zhǔn)偏差，從而得到目標(biāo)地?zé)崽锏臏囟?、體積、能量平均值和標(biāo)準(zhǔn)差。該方法的優(yōu)點是可給出待評價資源量的極值范圍、取值置信度和概率分布，提供預(yù)測值范圍和最樂觀值、最保守值等數(shù)據(jù)(Rybach et al.，1986)。

在使用熱儲法進(jìn)行地?zé)豳Y源評估時，一般將參數(shù)設(shè)為定值進(jìn)行計算。不同勘查階段、不同地?zé)崽铽@得的熱儲參數(shù)值通常為一個范圍，即在一定范圍內(nèi)變換的隨機變量。在進(jìn)行資源量評估時，熱儲的厚度、面積、溫度、孔隙度等主要參數(shù)均為范圍值和可能值，其精度決定了資源量評估結(jié)果的精度。蒙特卡羅法主要基于熱儲法(體積法)設(shè)定隨機變量和預(yù)測結(jié)果的性質(zhì)、范圍，以概率及統(tǒng)計學(xué)的理論方法為基礎(chǔ)，根據(jù)待求隨機問題的變化規(guī)律和物理現(xiàn)象自身的統(tǒng)計規(guī)律構(gòu)造概率模型，然后通過統(tǒng)計實驗使某些統(tǒng)計參量恰為待求問題的解(李同彪，2015)。

蒙特卡羅法在模擬計算地?zé)豳Y源量(預(yù)測發(fā)電量)時，一般通過設(shè)定隨機過程反復(fù)生成時間序列，設(shè)定參數(shù)估計量和統(tǒng)計量(統(tǒng)計結(jié)果)，同時給出結(jié)果的概率分布模式(函數(shù))。該方法可利用較多的熱儲資料，較好地解決計算結(jié)果的可信度問題(Nathenson，1978)。

模擬前應(yīng)設(shè)計模型(模擬公式)，設(shè)定公式中的隨機變量(變量)，提出問題(預(yù)測量)，設(shè)定變量和預(yù)測量的概率分布模型(分布函數(shù))；然后給出各隨機變量的值域和最可能取值，設(shè)計某種規(guī)則(模擬公式)，根據(jù)模型中各隨機變量的分布類型，通過大數(shù)據(jù)量的模擬試驗(隨機抽樣試驗)對隨機抽樣數(shù)據(jù)進(jìn)行數(shù)學(xué)計算、統(tǒng)計分析模擬試驗，使預(yù)測量的解對應(yīng)于模型中隨機變量的某些特征(如概率、均值等)，得出的預(yù)測量通過概率分布進(jìn)行顯示(顏英軍，2009；李同彪，2015)。

試驗結(jié)果包括最小值、最大值、數(shù)學(xué)期望值和標(biāo)準(zhǔn)偏差等，生成概率分布曲線圖以定義結(jié)果的類型(李同彪，2015)。

2.2.2 應(yīng)用實例 目標(biāo)地?zé)崽餅槁∑鹕降貙α餍偷責(zé)豳Y源，屬于變質(zhì)-火山巖基巖裂隙型熱儲。預(yù)測熱儲埋深在800～2 000 m之間，可視為圈閉型熱儲。基于熱儲法，將公式(1)變換以適于蒙特卡羅法的應(yīng)用(顏英軍，2009；韓征等，2015)。

Q儲=Q容=XAHφ

(3)

式(3)中：Q儲為地?zé)崃黧w存儲量，m3；Q容為地?zé)崃黧w容積儲量，m3；A為熱儲面積，m2；H為熱儲層厚度，m；φ為熱儲層孔隙度；X為可采系數(shù)，受熱儲層巖性、厚度、孔隙度等控制，淺層水的開采系數(shù)X取值一般在1%～2%之間，深層水X值稍大于淺層水(邱楠生等，2004)。

利用蒙特卡羅法模擬，設(shè)定95%置信度條件，給出100%、95%、75%、50%概率條件下的預(yù)測資源量(可發(fā)電量)。熱儲法模擬計算需要的參數(shù)(條件)有熱儲面積和厚度、熱儲參數(shù)、滲透率，以及基礎(chǔ)巖石學(xué)、熱力學(xué)、水力學(xué)參數(shù)等。

(1) 模擬公式及變量預(yù)測量定義。為研究地?zé)崽锏責(zé)豳Y源可利用量(發(fā)電量)，基于公式(3)設(shè)計蒙特卡羅法模擬計算的預(yù)設(shè)變量和預(yù)測變量，計算公式為：

(4)

式(4)中：Q為實際可發(fā)電量，MW；k為綜合利用系數(shù)，根據(jù)地?zé)岵墒章屎桶l(fā)電效率綜合而得；A為預(yù)測的熱儲面積，m2；H有為熱儲有效厚度，m；C綜為體積法計算的綜合熱容比，kJ/(kg·℃)；Tpg為熱儲溫度，℃；Tcut-off為根據(jù)地?zé)岚l(fā)電利用確定的地?zé)崃黧w下限溫度，℃，取120 ℃(印度尼西亞能源和礦產(chǎn)資源部，2011)；(30×7 000×3 600)為地?zé)岚l(fā)電時間(30年發(fā)電期，每年利用時長7 000 h)，s;×1 000表示將發(fā)電量結(jié)果換算成MW。

表1 B地?zé)崽锩商乜_法模擬計算參數(shù)值

H有=αH

(5)

式(5)中:α為儲厚比，參照基巖裂隙型熱儲的經(jīng)驗值，介于0.02～0.40之間；H為熱儲層所在地層總厚度，m。

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

(6)

式(6)中：ρw為T溫度時水的密度，kg/m3，估算熱儲溫度在100～200 ℃之間，取定值920；ρr為T溫度時巖石的密度，kg/m3，熱儲地層為變質(zhì)巖和侵入巖，取定值2 500；φ為孔隙度，根據(jù)鉆井實測值確定取值范圍為0.04～0.19；Cw為水熱容，kJ/(kg·℃)，取定值4.18；Cr為巖石熱容，kJ/(kg·℃)，采用定值0.879。

模擬時定義了資源量計算主要參數(shù)的函數(shù)類型(三角函數(shù)、正太函數(shù)、對數(shù)函數(shù)等)，確定合理的取值(最大值、最小值、平均值、眾數(shù)值等)范圍，參考地?zé)岬刭|(zhì)條件對模擬結(jié)果進(jìn)行分析，保證結(jié)果具有較好的可靠性。

發(fā)電量Q根據(jù)公式(4)、(5)、(6)模擬計算而得，式中采用的變量、預(yù)測量及范圍見表1。

(2) 模擬結(jié)果。使用Oracle Crystal Ball軟件,設(shè)置300 001次隨機模擬次數(shù)，計算B地?zé)崽锩商乜_法的預(yù)測結(jié)果，設(shè)置為95%的置信度區(qū)間，預(yù)計30年發(fā)電期，每年利用7 000 h，發(fā)電量單位為MW。具體結(jié)果如下。

① 100%概率條件下的預(yù)測值。圖3顯示了一種常用集中分布函數(shù)，屬于正態(tài)分布類型(正態(tài)函數(shù))，結(jié)果為對稱分布。最可能發(fā)生值是圖中顯示的平均值59.14 MW，標(biāo)準(zhǔn)偏差為30.28。圖3顯示，接近平均值的數(shù)發(fā)生幾率較高。

圖3 100%概率條件下B地?zé)崽镱A(yù)測發(fā)電量Fig. 3 Predicted electric generating capacity of the B geothermal field under the probability of 100%

根據(jù)切比雪夫定理(王學(xué)仁，1982)，任何樣本中，落在平均值的k個標(biāo)準(zhǔn)偏差范圍內(nèi)的觀測值例(幾率)至少為1-1/k2，其中k>1。當(dāng)k=3時，幾率>8/9，因此以平均值±3倍標(biāo)準(zhǔn)偏差為預(yù)測值的最大值，取值范圍為0～149.98 MW。

② 99%概率條件下的預(yù)測值。圖4顯示了正態(tài)分布函數(shù)。最可能發(fā)生值是平均值59.14 MW，在該概率下預(yù)測的最小值為8.92 MW，最大值為157.56 MW。

圖4 99%概率條件下B地?zé)崽镱A(yù)測發(fā)電量Fig. 4 Predicted electric generating capacity of the B geothermal field under the probability of 99%

③ 95%概率條件下的預(yù)測值。圖5顯示了一種正態(tài)分布函數(shù)。最可能發(fā)生值是平均值59.14 MW，該概率下預(yù)測的最小值為14.25 MW，最大值為129.55 MW。

圖5 95%概率條件下B地?zé)崽镱A(yù)測發(fā)電量Fig. 5 Predicted electric generating capacity of the B geothermal field under the probability of 95%

④ 75%概率條件下的預(yù)測值。圖6顯示了正態(tài)分布函數(shù)。最可能發(fā)生值是平均值59.14 MW，在該概率下預(yù)測的最小值為26.18 MW，最大值為95.26 MW。

圖6 75%概率條件下B地?zé)崽镱A(yù)測發(fā)電量Fig. 6 Predicted electric generating capacity of the B geothermal field under the probability of 75%

⑤ 50%概率條件下預(yù)測值。圖7顯示了正態(tài)分布函數(shù)。最可能發(fā)生值是平均值59.14 MW，該概率下預(yù)測的最小值為36.3 MW，最大值為77.2 MW。

圖7 50%概率條件下B地?zé)崽镱A(yù)測發(fā)電量Fig. 7 Predicted electric generating capacity of the B geothermal field under the probability of 50%

⑥ 隨機變量權(quán)重分析。B地?zé)崽镱A(yù)測模型敏感度圖解(圖8)顯示，熱儲面積參數(shù)V2占比55.5%，熱儲溫度參數(shù)V6占比37.6%，綜合利用系數(shù)參數(shù)V8占比3.4%，有效熱儲厚度參數(shù)V3占比3.3%，綜合熱導(dǎo)系數(shù)參數(shù)V7占比0.2%。

敏感度圖顯示了模擬公式中隨機變量參數(shù)對結(jié)果的影響程度。圖8表明，熱儲面積、熱儲溫度、利用系數(shù)、熱儲厚度、綜合導(dǎo)熱系數(shù)5項參數(shù)為模擬結(jié)果(發(fā)電量)的主要影響因子，且對結(jié)果的影響程度依次降低。

圖8 B地?zé)崽镱A(yù)測模型敏感度圖解Fig. 8 Sensitive chart of predicted model of the B geothermal field

3 討 論

3.1 預(yù)測發(fā)電量最小值

圖3—圖7為蒙特卡羅法在不同概率下的結(jié)果和函數(shù)分布類型，根據(jù)正態(tài)函數(shù)性質(zhì)，100%概率模式下的最小值>0，但未給出實際值。

根據(jù)地表熱釋放量法計算的結(jié)果為834.86 kW(印度尼西亞能源和礦產(chǎn)資源部，2011)，根據(jù)世界范圍內(nèi)的經(jīng)驗值(中國科學(xué)院地質(zhì)研究所地?zé)峤M，1978)，地面釋放量一般比實際值至少差1～3個數(shù)量級(李同彪，2015)，以10倍對其修正后得出的值(8.348 6 MW)作為該方法可利用資源量的最小下限值。

運用體積法，設(shè)計發(fā)電期30年、年利用時長為7 000 h、SiO2的熱儲溫度為161 ℃、回收溫度為120 ℃、孔隙度為10%，估算資源量為27.00 MW。

兩者對比，體積法更接近且小于實際預(yù)測的發(fā)電量，因此以27.00 MW為預(yù)測發(fā)電量下限值。

3.2 蒙特卡羅法預(yù)測值

蒙特卡羅法可用于地?zé)衢_發(fā)各階段資源量的預(yù)測(王心義等，2002；韓征等，2015)，對有多年動態(tài)監(jiān)測資料的地?zé)崽锏馁Y源量估算更準(zhǔn)確(顏英軍，2009)，在估算資源量時可給出置信區(qū)間(李同彪，2015)。研究區(qū)尚無地?zé)徙@孔，有效的地?zé)豳Y源計算方法為體積法，但其計算參數(shù)均為定值，無法給出置信區(qū)間，因此采用已有資料和分析數(shù)據(jù)時均選用保守參數(shù)值來計算，可認(rèn)為在目前地?zé)嵴{(diào)查階段59.14 MW的地?zé)岚l(fā)電量可信度高。隨著勘查工作的深入，其資源量呈增加趨勢。

4 結(jié) 論

(1) 此次地?zé)豳Y源量預(yù)測僅針對1 km以內(nèi)的淺層熱儲，熱儲面積僅為B溫泉所在的小片區(qū)域。

(2) 結(jié)合體積法固定參數(shù)計算和模擬試驗結(jié)果分析，認(rèn)為資源量計算的最可能值在置信度為95%時為59.14 MW，可能范圍為27.00～149.98 MW，遠(yuǎn)景發(fā)電量為115.85 MW。

(3) 影響發(fā)電量預(yù)測值的主要參數(shù)隨機率降低依次為預(yù)測熱儲面積、熱儲溫度、利用系數(shù)、熱儲厚度、綜合導(dǎo)熱系數(shù)。

(4) 根據(jù)大地構(gòu)造環(huán)境、構(gòu)造巖漿演化及鄰區(qū)鉆孔資料分析，研究區(qū)深部存在火山巖漿熱源的裂隙型熱儲，特別是B溫泉區(qū)深部存在>200 ℃的基巖裂隙型熱儲，具有有利的變質(zhì)巖熱儲蓋層和由斷裂構(gòu)造溝通的地?zé)崃黧w循環(huán)、補給系統(tǒng)，在深度>1 km的深部有第二熱儲，具有較大的地?zé)豳Y源開發(fā)潛力。