數值模擬在地熱資源評價及優化開采方案的應用研究

魏永霞 ，唐仲華 ，左霖

(1安徽省地質環境監測總站， 安徽合肥 230001 ；2中國地質大學， 湖北武漢 430074 )

0 引言

昭關溫泉蘊藏于安徽省巢湖市含山縣昭關鎮轄區，即當年“伍子胥過昭關，一夜愁白了頭”的歷史名勝之地，距含山縣城約10.8km。含山縣未來發展定位：長三角休閑度假旅游目的地，而溫泉的加入則會大大提升旅游的品質，因此，摸清昭關地熱資源量及開采現狀、未來開采潛力成為當務之急。在常規勘查的基礎上，采用數值模擬則能較好地解決地熱田合理開發利用的問題。

1 地熱田地質背景

昭關地處亞熱帶濕潤季風氣候區，地形東南高西北低，最高為東部的雙尖山256.6m，其次為南部虎山196.2m，最低處為大馬村一帶8.5m。

區內地層有震旦系-寒武系、志留系及白堊系、第四系：震旦系僅出露燈影組，厚度＞164m，巖性為白云巖；寒武系發育齊全，厚度30至190m不等，巖性亦為白云巖；志留系為高家邊組、墳頭組，厚度均＞300m，巖性以砂巖、泥巖、頁巖為主；白堊系包括浦口組、赤山組，厚度＞50m，巖性主要為礫巖、細砂巖、粉砂巖、粉砂質泥巖；第四系厚度一般10～20m，巖性為粉質黏土。巖漿巖分布在地熱田東側，呈橢圓狀，南北長軸約450m，東西短軸約350m，巖體為淺火山巖。

構造為香泉推覆體，以震旦系—奧陶系所組成的大型平臥褶皺線性緊閉倒轉（同斜）褶皺及一系列疊瓦狀逆沖斷層發育，北東向推覆斷層F0、F3斷層、F4斷層及北西向F11斷層、F12斷層發育。地熱即出露于香泉推覆體的后緣的傾沒端，如圖1。

圖1 地熱田地質圖Figure 1. Geological map of the geothermal field

2 地熱田概念模型

地熱田平面展布：F4斷層西北側志留系砂頁巖，對地熱流體起到阻隔作用，確定為地熱隔水邊界；F11斷層東北側和F12斷層西南側無地熱排泄顯示，阻隔了地熱的排泄；東南部F0斷層為地熱補給通道。地熱田北東-南西長1.58km，北西-南東寬1.12km，面積1.241km2（圖2）。

圖2 地熱田熱儲概念模型Figure 2. Conceptual model of geothermal field reservoir

地熱田垂向展布：北部淺，一般12～160m；南部深，由西南至東北：上界面600～870～1170m，下界面一般在850～1300～1600m，為一向南東傾斜的斷面似“亞玲”狀的結構體，并且西南部為埋藏型。熱儲厚度一般在260～360m之間，排泄區碳酸鹽巖巖溶裂隙及破碎段地層總厚約180m，地熱田平均厚度約310m。

綜上所述，F0為地熱的補給邊界，F3為地熱的導熱儲熱邊界，F11和F12為弱透熱邊界，F4為地熱的隔熱邊界；玄武巖體為阻熱邊界，松散層和巖溶發育微弱的碳酸鹽巖為地熱的隔熱邊界。

圖3 地熱田三維立體圖Figure 3. Three-dimensional map of the geothermal field

3 地熱田模型的建立

3.1 地熱田三維地質模型建立

三維地質模型主要基于TINs生成的。根據地層巖性、斷層的展布及地面高程點等信息，通過Mapgis、Surfer、GMS軟件生成的一系列TINs建立了三維地質結構模型（如圖3）。

3.2 地熱田數學模型的建立

3.2.1 地下水流動數學模型

根據達西定律和水均衡原理建立了地下水流動的微分方程為：

初始條件為：

邊界條件為：

式中：Ω-滲流區域；H-地下水系統的水頭（m）；Kxx、Kyy、Kzz-含水層在滲透主方向上的滲透系數（m/d）；ε-含水層的源匯項（m/d）；H0-初始流場（m）；μ-比儲水系數（1/m）；Γ2—第二類邊界（已知流量的邊界）；Kn-第二類邊界法向方向的滲透系數（m/d），n-邊界面的外法線方向；q-第二類邊界Γ2上單位面積上的流量（m/d），流入為正，流出為負。

3.2.2 熱運移數學模型

對于中低溫地下熱水系統，地下熱水主要以液態的形式存在，因此在多孔介質中地下熱水熱量運移的模式主要分為三種：一是通過水的對流作用傳輸熱量；二是熱量通過多孔介質傳遞；三是通過水的熱機械彌散作用運輸熱量（薛禹群等，1990；薛禹群等，1994）。根據地下熱水的運移模式和Bear（1972）提出的能量守恒原理，得到熱運移方程為：

式中：ρw、ρc-地下水和巖石的密度（kg/m3）；Cw、Cc-地下水和巖石的比熱（J /kg·℃）；φ-有效孔隙度；T-地熱系統的溫度（℃）；qi-地下水滲流速度（m/d）；λi,j-熱傳導系數；Q-源匯項；ρQ、CQ-源匯項的密度和比熱；Tq-抽/注水井水的溫度。

圖4 地熱田三維剖分網格圖Figure 4. Three-dimensional grid diagram of the geothermal field

4 地熱田數值模擬

4.1 地熱田地下水流數值模擬

平面上，該地熱田東北部為敞開的補給邊界，F0、F3斷層概化為隔水邊界，垂向上深度為1600m，上部、底部均為零流量邊界。采用矩形網格對地熱田進行剖分，上下劃分8層，平面上剖分成100*100個網格，共有80000個單元，91809個節點，其中有效單元31688個，每層3961個有效單元。數值模型剖分網格如圖4。

根據地熱田區地熱井、溫泉水位監測資料，確定了模擬期及初始水位；根據自流情況及同位素分析，確定了模型的源匯項；根據鉆孔資料及抽水試驗成果，得到含水層基本水文地質參數，并在模型校正和檢驗中，不斷地調整，最終得出接近實際的參數值。整個模擬期內地下水總補給量為5.72×104m3/a，總排泄量為6.28×104m3/a，補排差為-0.56×104m3/a，為負均衡。故在地熱田區內應禁止人為開采。

4.2 地熱田熱對流-傳導數值模擬

熱運移模型的范圍、邊界條件及空間網格剖分與時間網格剖分均與水流模型一致。其中隔水邊界即為零熱流邊界，上部邊界為恒溫邊界，溫度為18.3℃，底部邊界根據研究區大地熱流和巖性計算熱流，平均為153.26 mW/m2。巖石的熱物性參數采用經驗值。根據監測孔溫度變化資料對模型進行檢驗校正，使模擬溫度接近實測溫度，如圖5和圖6。

圖5 模擬溫度分布圖Figure 5. Simulated temperature profile

圖6 實測溫度分布圖Figure 6. Distribution of measured temperatures

5 地熱田地熱資源量

地熱勘查程度、成熱條件、資源類型、已有資料豐富程度等不同，資源量的計算方法亦不相同，根據地熱田的特點，本文采用熱儲法和基于數值模擬的水熱均衡法兩種方法進行對比計算。

5.1 熱儲法計算的地熱資源量

根據昭關地熱田熱儲模型，地熱田是主要受深大活動斷裂控制的褶皺斷裂復合型熱儲，地熱資源計算區劃分為以層狀熱儲為主的層狀熱儲區和以斷裂帶為主的帶狀熱儲區，將平均熱儲層底板深度到熱礦水最大循環深度作為斷裂帶帶狀熱儲區的計算厚度，其橫向影響帶寬度作為熱儲寬度，根據各區熱儲特征研究確定所對應的熱儲參數進行資源量的計算。計算結果，地熱田熱儲平均熱容量253.0×104J/m3·℃，蘊藏的總資源量為4.62×1016J，其中巖石蘊藏的熱能為4.32×1016J，熱液蘊藏的熱能為3×1015J。地熱水儲存量9.74×106m3，地熱水儲存量換算成的熱量6.21×1015J，可采資源量8.32×1015J，計算結果如表1所示。

表1 地熱資源量計算結果Table 1. Calculation results of geothermal resources

表2 地熱資源量模擬成果表Table 2. Simulation results of geothermal resources

5. 2 基于數值模擬計算的地熱資源量

根據地熱水溫度將計算區劃分為大于60℃地區，50≤T＜60℃的地區， 25≤T＜50℃的地區。數值模擬計算結果見表2，地熱田熱儲中蘊藏的總資源量為5.67×1016J，其中巖石蘊藏的熱能為5.42×1016J，熱液蘊藏的熱能為2.53×1015J。地熱水儲存量15.18×106m3。

兩種方法計算結果對比可知，數值模擬結果比熱儲法的略大一些，主要因為熱儲法計算參數采用的是整個區域上的平均值，數值模擬時采用的是各個單元上的值，相對來說精確要高一些。

6 地熱田開采現狀

地熱田區地熱流體主要用于熱帶魚養殖，開采的地熱井及溫泉主要分布在玄武巖體的西南，為HS04、HS06、HS08及HS12井4開采井和ZJ06（大湖潭）、ZJ07（小湖潭）2個泉群，共6處，地熱的開采主要采用自流直排的方式，開采量與井（泉）自流量基本一致（圖7），2015年-2016年7月開采量為2958～5238m3/d，平均為4929.22 m3/d。

圖7 開采量與水位動態關系Figure 7. Dynamic relation between yield and water level

7 優化開采方案

通過檢驗校正后的模型，對地熱田的三維非穩定流進行預測模擬，以確定合理的開采井和開采量，為昭關地熱資源的可持續利用提供依據。

7.1 方案一

為使模型的開采方案更符合實際，將模型開采方案對模擬期2016年7月—2021年7月5年中豐水年、平水年、枯水年分三種不同組合進行模擬：

組合一：2年豐水年、3年平水年（簡稱“2豐3平”），側向補給量為現狀條件下增加20%，預測溫度場的變化情況如圖8、圖9。觀測點ob1水位前2年豐水年側向流量補給量大水位出現少許上升，隨著后3年平水年側向流量補給恢復正常，水位回落，但略高于初始水位；觀測點HS08溫度在29℃上下波動。

圖8 “2豐3平”條件下觀測點ob1水位變化曲線圖Figure 8. Curve of water level change at the observation point ob1 under the condition of "2 wet and 3 normal years"

圖9 “2豐3平”條件下觀測點HS08溫度變化曲線圖Figure 9. Temperature change curve of the observation point HS08 under the condition of “2 wet and 3 normal years”

組合二：1年為豐水年，1年為枯水年，3年為平水年（簡稱“1豐1枯3平”），且側向流量補給保持現狀條件不變，預測溫度場的變化情況如圖10、圖11所示。觀測點ob1前期豐水年側向流量補給量大水位出現少許上升，隨著枯水年側向流量補給減少，水位下降，后期平水年側向流量補給回升，水位出現緩慢上升趨勢，最終基本與初始水位持平；觀測點HS08溫度在29℃上下波動。

圖10 “1豐1枯3平”條件下觀測點ob1水位變化曲線圖Figure 10. Water level variation curve at the observation point ob1 under the condition of "1 wet, 1 dry and 3 normal years"

圖11 “1豐1枯3平”條件下觀測點HS08溫度變化曲線圖Figure 11. Temperature change curve of the observation point HS08 under the condition of "1 wet, 1 dry and 3 normal"

組合三：2年為枯水年，3年為平水年（簡稱“2枯3平”），且側向流量補給為現狀條件下減少20%，預測溫度場的變化情況如圖12、圖13所示。觀測點ob1水位前期下降后期略有上升，最終略高于初始水位；觀測點HS08溫度在29℃上下波動。

7.2 方案二

圖12 “2枯3平”條件下觀測點ob1水位變化曲線圖Figure 12. Curve of water level at the observation point ob1 under the condition of "2 dry and 3 normal years"

圖13 “2枯3平”條件下觀測點HS08溫度變化曲線圖Figure 13. Temperature change curve of the observation point HS08 under the condition of "2 dry and 3 normal years"

由于現狀開采條件下平水期、枯水期水位均出現下降趨勢，為了使開采井和泉保持可持續開采，因此需要減少現狀開采量使地下水水位不出現下降情況，經過不斷地調試開采量的大小得出正常年份下現狀開采量減少15 %的情況下，地下水水位基本保持不變，溫度場也無異常。(如圖14、圖15)。

圖14 觀測點HS04水位變化曲線圖Figure 14. Curve of water level change at the observation point HS04

圖15 觀測點HS08溫度變化曲線圖Figure 15. Temperature change curve of the observation point HS08

對比方案預測結果，綜合分析得出方案二更為合理，5年中地下水水位沒有出現下降趨勢，溫度變化比較穩定。

8 結論

通過建立地熱系統水-熱數值模擬模型，計算出地熱田熱儲中蘊藏的總資源量為5.6728×1016J，其中巖石蘊藏的熱能為5.4192×1016J，熱液蘊藏的熱能為2.5302×1015J。地熱水儲存量15.18×106m3。

擬定了4中不同開采方案，并利用數值模型進行了預測評價，得出優化方案即現狀開采量減少15%，地下水水位沒有出現下降趨勢，溫度場也無異常。