数值模拟在地热资源评价及优化开采方案的应用研究

魏永霞 ，唐仲华 ，左霖

(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%，地下水水位没有出现下降趋势，温度场也无异常。