數(shù)值模擬在地?zé)豳Y源評價及優(yōu)化開采方案的應(yīng)用研究

魏永霞 ，唐仲華 ，左霖

(1安徽省地質(zhì)環(huán)境監(jiān)測總站， 安徽合肥 230001 ；2中國地質(zhì)大學(xué)， 湖北武漢 430074 )

0 引言

昭關(guān)溫泉蘊藏于安徽省巢湖市含山縣昭關(guān)鎮(zhèn)轄區(qū)，即當(dāng)年“伍子胥過昭關(guān)，一夜愁白了頭”的歷史名勝之地，距含山縣城約10.8km。含山縣未來發(fā)展定位：長三角休閑度假旅游目的地，而溫泉的加入則會大大提升旅游的品質(zhì)，因此，摸清昭關(guān)地?zé)豳Y源量及開采現(xiàn)狀、未來開采潛力成為當(dāng)務(wù)之急。在常規(guī)勘查的基礎(chǔ)上，采用數(shù)值模擬則能較好地解決地?zé)崽锖侠黹_發(fā)利用的問題。

1 地?zé)崽锏刭|(zhì)背景

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

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

構(gòu)造為香泉推覆體，以震旦系—奧陶系所組成的大型平臥褶皺線性緊閉倒轉(zhuǎn)（同斜）褶皺及一系列疊瓦狀逆沖斷層發(fā)育，北東向推覆斷層F0、F3斷層、F4斷層及北西向F11斷層、F12斷層發(fā)育。地?zé)峒闯雎队谙闳聘搀w的后緣的傾沒端，如圖1。

圖1 地?zé)崽锏刭|(zhì)圖Figure 1. Geological map of the geothermal field

2 地?zé)崽锔拍钅Ｐ?/h2>

地?zé)崽锲矫嬲共迹篎4斷層西北側(cè)志留系砂頁巖，對地?zé)崃黧w起到阻隔作用，確定為地?zé)岣羲吔纾籉11斷層?xùn)|北側(cè)和F12斷層西南側(cè)無地?zé)崤判癸@示，阻隔了地?zé)岬呐判梗粬|南部F0斷層為地?zé)嵫a給通道。地?zé)崽锉睎|-南西長1.58km，北西-南東寬1.12km，面積1.241km2（圖2）。

圖2 地?zé)崽餆醿Ω拍钅Ｐ虵igure 2. Conceptual model of geothermal field reservoir

地?zé)崽锎瓜蛘共迹罕辈繙\，一般12～160m；南部深，由西南至東北：上界面600～870～1170m，下界面一般在850～1300～1600m，為一向南東傾斜的斷面似“亞玲”狀的結(jié)構(gòu)體，并且西南部為埋藏型。熱儲厚度一般在260～360m之間，排泄區(qū)碳酸鹽巖巖溶裂隙及破碎段地層總厚約180m，地?zé)崽锲骄穸燃s310m。

綜上所述，F(xiàn)0為地?zé)岬难a給邊界，F(xiàn)3為地?zé)岬膶?dǎo)熱儲熱邊界，F(xiàn)11和F12為弱透熱邊界，F(xiàn)4為地?zé)岬母魺徇吔纾恍鋷r體為阻熱邊界，松散層和巖溶發(fā)育微弱的碳酸鹽巖為地?zé)岬母魺徇吔纭?/p>

圖3 地?zé)崽锶S立體圖Figure 3. Three-dimensional map of the geothermal field

3 地?zé)崽锬Ｐ偷慕?/h2>

3.1 地?zé)崽锶S地質(zhì)模型建立

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

3.2 地?zé)崽飻?shù)學(xué)模型的建立

3.2.1 地下水流動數(shù)學(xué)模型

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

初始條件為：

邊界條件為：

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

3.2.2 熱運移數(shù)學(xué)模型

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

式中：ρw、ρc-地下水和巖石的密度（kg/m3）；Cw、Cc-地下水和巖石的比熱（J /kg·℃）；φ-有效孔隙度；T-地?zé)嵯到y(tǒng)的溫度（℃）；qi-地下水滲流速度（m/d）；λi,j-熱傳導(dǎo)系數(shù)；Q-源匯項；ρQ、CQ-源匯項的密度和比熱；Tq-抽/注水井水的溫度。

圖4 地?zé)崽锶S剖分網(wǎng)格圖Figure 4. Three-dimensional grid diagram of the geothermal field

4 地?zé)崽飻?shù)值模擬

4.1 地?zé)崽锏叵滤鲾?shù)值模擬

平面上，該地?zé)崽飽|北部為敞開的補給邊界，F(xiàn)0、F3斷層概化為隔水邊界，垂向上深度為1600m，上部、底部均為零流量邊界。采用矩形網(wǎng)格對地?zé)崽镞M行剖分，上下劃分8層，平面上剖分成100*100個網(wǎng)格，共有80000個單元，91809個節(jié)點，其中有效單元31688個，每層3961個有效單元。數(shù)值模型剖分網(wǎng)格如圖4。

根據(jù)地?zé)崽飬^(qū)地?zé)峋?、溫泉水位監(jiān)測資料，確定了模擬期及初始水位；根據(jù)自流情況及同位素分析，確定了模型的源匯項；根據(jù)鉆孔資料及抽水試驗成果，得到含水層基本水文地質(zhì)參數(shù)，并在模型校正和檢驗中，不斷地調(diào)整，最終得出接近實際的參數(shù)值。整個模擬期內(nèi)地下水總補給量為5.72×104m3/a，總排泄量為6.28×104m3/a，補排差為-0.56×104m3/a，為負均衡。故在地?zé)崽飬^(qū)內(nèi)應(yīng)禁止人為開采。

4.2 地?zé)崽餆釋α?傳導(dǎo)數(shù)值模擬

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

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

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

5 地?zé)崽锏責(zé)豳Y源量

地?zé)峥辈槌潭?、成熱條件、資源類型、已有資料豐富程度等不同，資源量的計算方法亦不相同，根據(jù)地?zé)崽锏奶攸c，本文采用熱儲法和基于數(shù)值模擬的水熱均衡法兩種方法進行對比計算。

5.1 熱儲法計算的地?zé)豳Y源量

根據(jù)昭關(guān)地?zé)崽餆醿δＰ停責(zé)崽锸侵饕苌畲蠡顒訑嗔芽刂频鸟薨檾嗔褟?fù)合型熱儲，地?zé)豳Y源計算區(qū)劃分為以層狀熱儲為主的層狀熱儲區(qū)和以斷裂帶為主的帶狀熱儲區(qū)，將平均熱儲層底板深度到熱礦水最大循環(huán)深度作為斷裂帶帶狀熱儲區(qū)的計算厚度，其橫向影響帶寬度作為熱儲寬度，根據(jù)各區(qū)熱儲特征研究確定所對應(yīng)的熱儲參數(shù)進行資源量的計算。計算結(jié)果，地?zé)崽餆醿ζ骄鶡崛萘?53.0×104J/m3·℃，蘊藏的總資源量為4.62×1016J，其中巖石蘊藏的熱能為4.32×1016J，熱液蘊藏的熱能為3×1015J。地?zé)崴畠Υ媪?.74×106m3，地?zé)崴畠Υ媪繐Q算成的熱量6.21×1015J，可采資源量8.32×1015J，計算結(jié)果如表1所示。

表1 地?zé)豳Y源量計算結(jié)果Table 1. Calculation results of geothermal resources

表2 地?zé)豳Y源量模擬成果表Table 2. Simulation results of geothermal resources

5. 2 基于數(shù)值模擬計算的地?zé)豳Y源量

根據(jù)地?zé)崴疁囟葘⒂嬎銋^(qū)劃分為大于60℃地區(qū)，50≤T＜60℃的地區(qū)， 25≤T＜50℃的地區(qū)。數(shù)值模擬計算結(jié)果見表2，地?zé)崽餆醿χ刑N藏的總資源量為5.67×1016J，其中巖石蘊藏的熱能為5.42×1016J，熱液蘊藏的熱能為2.53×1015J。地?zé)崴畠Υ媪?5.18×106m3。

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

6 地?zé)崽镩_采現(xiàn)狀

地?zé)崽飬^(qū)地?zé)崃黧w主要用于熱帶魚養(yǎng)殖，開采的地?zé)峋皽厝饕植荚谛鋷r體的西南，為HS04、HS06、HS08及HS12井4開采井和ZJ06（大湖潭）、ZJ07（小湖潭）2個泉群，共6處，地?zé)岬拈_采主要采用自流直排的方式，開采量與井（泉）自流量基本一致（圖7），2015年-2016年7月開采量為2958～5238m3/d，平均為4929.22 m3/d。

圖7 開采量與水位動態(tài)關(guān)系Figure 7. Dynamic relation between yield and water level

7 優(yōu)化開采方案

通過檢驗校正后的模型，對地?zé)崽锏娜S非穩(wěn)定流進行預(yù)測模擬，以確定合理的開采井和開采量，為昭關(guān)地?zé)豳Y源的可持續(xù)利用提供依據(jù)。

7.1 方案一

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

組合一：2年豐水年、3年平水年（簡稱“2豐3平”），側(cè)向補給量為現(xiàn)狀條件下增加20%，預(yù)測溫度場的變化情況如圖8、圖9。觀測點ob1水位前2年豐水年側(cè)向流量補給量大水位出現(xiàn)少許上升，隨著后3年平水年側(cè)向流量補給恢復(fù)正常，水位回落，但略高于初始水位；觀測點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平”），且側(cè)向流量補給保持現(xiàn)狀條件不變，預(yù)測溫度場的變化情況如圖10、圖11所示。觀測點ob1前期豐水年側(cè)向流量補給量大水位出現(xiàn)少許上升，隨著枯水年側(cè)向流量補給減少，水位下降，后期平水年側(cè)向流量補給回升，水位出現(xiàn)緩慢上升趨勢，最終基本與初始水位持平；觀測點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平”），且側(cè)向流量補給為現(xiàn)狀條件下減少20%，預(yù)測溫度場的變化情況如圖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"

由于現(xiàn)狀開采條件下平水期、枯水期水位均出現(xiàn)下降趨勢，為了使開采井和泉保持可持續(xù)開采，因此需要減少現(xiàn)狀開采量使地下水水位不出現(xiàn)下降情況，經(jīng)過不斷地調(diào)試開采量的大小得出正常年份下現(xiàn)狀開采量減少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

對比方案預(yù)測結(jié)果，綜合分析得出方案二更為合理，5年中地下水水位沒有出現(xiàn)下降趨勢，溫度變化比較穩(wěn)定。

8 結(jié)論

通過建立地?zé)嵯到y(tǒng)水-熱數(shù)值模擬模型，計算出地?zé)崽餆醿χ刑N藏的總資源量為5.6728×1016J，其中巖石蘊藏的熱能為5.4192×1016J，熱液蘊藏的熱能為2.5302×1015J。地?zé)崴畠Υ媪?5.18×106m3。

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