通州区地热资源优化开采模式动态研究

张进平，王煜曦，刘桂宏，何铁柱，袁利娟，徐浩然，孔祥军

(1. 北京市地热调查研究所，北京 102218；2. 清华大学土木工程系，北京 100084)

热储工程数值模拟技术已在地热资源评估、开采回灌方案优化设计以及地热井开采寿命预测等方面发挥重要作用。在地热开发早期，由于已知数据较少，采用高度理想化、概念清晰、计算量相对较小的解析模型[1]或集中参数模型[2]对热储层特性演化特征进行粗略预估。随着计算机技术的发展，数值模拟可以精细刻画热储层非均质性及其热储物理、化学、热力学等在开采回灌过程中的演化规律，数值模拟已成为现在热储工程中多场耦合效应研究与评价的主要方法[3]。热储工程热-水-力-化(THMC)多场耦合模型主要有两类：连续介质力学模型和离散裂隙网络模型[4-6]。基于连续介质力学模型的热储工程数值模拟方法发展很快，其中以有限元方法应用最为广泛。由于地热井(直径为0.1 m～0.2 m)和地热田(面积大于100 km2)之间的尺度相差巨大，群井系统数值模拟通常无法在日常工作计算机上完成，而城市尺度地热田群井系统高效模拟方法解决了这一难题[7]。目前热储模型大多是对井系统，李腾风等[8]考虑热源作用下非饱和土体水热耦合作用机制，基于格子Boltzmann方法，采用双分布函数分别描述温度场及水分场的演化过程，建立了相应的水热耦合模型。徐浩然等[9]提出了一种模拟工程尺度碳酸盐岩热储酸化压裂过程的数值方法，考虑碳酸盐岩热储酸化压裂过程中热、水、力、化四场之间的耦合作用。地热田含有数十口以上的地热井，故考虑群井效应对地热田可持续开发有着重要意义。通州区为北京市行政副中心，其地热勘探数据较为完整，为了确保其地热资源优势在本区内得到最大程度的发挥，构建科学高效的地热资源开发利用体系，文章结合了实际地层和断裂带构建了通州区地质模型，开展深层热储数值模拟，探寻北京市通州区深部地热资源的优化开采模式，为通州区规模化地热资源开发利用提供依据。

1 地质背景及热储模型建立

1.1 地质背景和地热井分布

野外实地调查及资料显示[10-11]，通州全区地势平坦而略有起伏，主要地貌类型为冲积平原，由永定河、潮白河、温榆河等冲洪积而成，可划分为河流阶地、河床及漫滩、砂质决口扇及砂丘、洼地、采土和砂石坑等5 种类型。

工作区基岩地质构造如图1，工作区新生界下伏地层从老到新包括中上元古界、古生界和中生界地层，主要发育夏垫断裂带、西集断裂、姚辛庄断裂等多条断裂带，其中夏垫断裂带对工作区钻井影响最大, 其余断裂规模相对较小，对地层的错动、切割等影响较小，对地层盖层条件、深部储层埋深、储层温度等影响相对较小。

图1 北京市通州区基岩地质构造图Fig. 1 Geological structure of bedrock in Beijing Tongzhou District

区内主要地层从新到老依次为新生界、石炭-二叠系、寒武-奥陶系、青白口系和蓟县系。新生界地层岩性以砂质粘土、粉细砂、泥岩、页岩为主，总厚度巨大，超过2000 m。新生界之下为古生界石炭-二叠系、奥陶系、寒武系，主要岩性为泥岩、页岩、砂岩、灰岩夹页岩。古生界之下为元古界地层，包括青白口系和蓟县系等；青白口系岩性以泥灰岩、砂岩、页岩为主；蓟县系包括铁岭组、洪水庄组、雾迷山组和杨庄组，铁岭组、雾迷山组和杨庄组岩性都以白云岩为主，洪水庄组岩性为页岩。根据地层岩性，工作区内主要以蓟县系雾迷山组热储为主。

根据以上资料建立通州区三维地质模型，区域面积 906 km2，深度4.5 km，地层分为 7 层，包含20 条断裂带，如图2(a)。地层由新到老依次为新生界(Kz)、石炭—二叠系(C-P)、寒武—奥陶系(Є-O)、青白口(Qn)、蓟县系铁岭组和洪水庄组(Jxt+Jxh)、蓟县系雾迷山组(Jxw)和太古界(Ar)如图2(b)。

图2 通州区三维地质模型示意图Fig. 2 Three-dimensional geological model of Tongzhou District

新生界全区范围内均有分布，通州西南部的地层厚度超过了4500 m。石炭—二叠系只在通州南部小部分区域存在，其余地区缺失。寒武—奥陶系在通州区的东部和南部分布较完整，北部大部分地层缺失，只在副中心附近小区域存在。青白口系、蓟县系铁岭组和洪水庄组在工作区中部及西南部有不同程度的缺失，根据之前的钻井资料，洪水庄组厚度约70 m 左右，为了简化模型，将铁岭组和洪水庄组合并成一层。蓟县系雾迷山组为通州区主要开采的热储层，在全区范围内均有分布，在南部及东部有部分地区缺失。太古界为基岩，全区范围内均有分布，由于钻井并未揭穿雾迷山组，东南部地区的太古界地层埋深可能超过4500 m，地层厚度分布及模型物理力学参数如表1 所示。

表1 地层分布及物理力学参数Table 1 Formation distribution and mechanical parameters

通州钻凿完成地热井 25 眼，其中开采井 22 眼，回灌井3 眼，主要集中在副中心区域，南部及西部处于未开发状态，没有地热井，地热井分布如图3。地热井成井深度、出水温度和出水量差异较大，现有数据统计：成井深度约642 m～3589 m，出水温度约28 ℃～92 ℃，出水量约603 m3/d～3073 m3/d。

图3 地热井分布图Fig. 3 Distribution of geothermal wells

1.2 多尺度模型的建立

COMSOL Multiphysics 是一款大型多物理场仿真软件，以有限元法为基础，通过求解偏微分方程、偏微分方程组来解决单物理场、多物理场问题，实现真实物理系统的模拟。多物理场问题本质上就是对偏微分方程组的求解，因此只要可以用偏微分方程组进行描述的物理现象，都可以用COMSOL Multiphysics 多物理场耦合来进行模拟分析。强大的多物理场耦合功能是该软件相比于其他有限元软件的主要优势所在。

1.2.1 地热井模型简化

李馨馨等[12]利用线单元模拟地热对井提出了裂隙岩体三维热流耦合的等效模拟方法，李敬元等[13]将地层岩石视为多孔介质，考虑了孔隙压力和渗流作用的影响，根据弹性理论和摩尔-库仑屈服条件以及岩石报伤后的软化特性，对井筒周围岩石进行了弹塑性力学分析，给出了井壁稳定条件。应用COMSOL 软件建立远场尺度热储三维数值模型，对于细长的地热井结构，其长度与半径之比可达104以上，是建模中挑战性的数值问题，Al-Khoury 等[14-15]和Saeid 等[16]提出了一种用一维线单元模拟地热井中传热过程的伪三维井筒模型。

井筒结构示意图如图4(a)所示，盖层以上部分放入套管，并采用水泥砂浆完井。本文进一步扩展了此模型，引入一维线单元对该三维井筒结构进行简化[17-18]，图4(a)为地热井筒结构示意图，图4(b)为简化的一维井筒模型，考虑沿井筒轴向的渗流传热过程，井筒内流体与盖层的热交换通过等效换热系数来近似考虑，储层部分放入滤水管，提供流体交换通道，其它主要假设条件包括：

图4 井筒结构和一维地热井模型Fig. 4 Well structure and one-dimensional model

1)储层处于完全饱和状态，且不考虑储层内气液相变过程；

2)盖层处于完全干燥状态，且不考虑不同储层之间的水力联系；

3)井筒内流体沿轴向流动，且同一深度流速沿径向处处相等；

4)考虑井筒内流体性质如密度、粘度、热传导系数、比热容等与温度的相关性，但不考虑井筒内气液相变过程。

1.2.2 模型网格划分与时间离散

兼顾计算精度与效率，地热井及周围区域的网格细化，最大单元尺寸2.5 m；热储层区域网格细化，最大单元尺寸1 km；其余盖层网格粗化，最大单元尺寸6 km，如图5 所示。模型总计包含299 959 个四面体单元，25 眼地热井由798 个一维线单元代表。

图5 模型网格剖分示意图Fig. 5 Mesh generation diagram of the model

地热系统只在每年的11 月15 日至来年的3 月15 日共4 个月的供暖季里处于开采或回灌运行状态，在非供暖季的8 个月里地热井处于停运状态，图6 为地热井运行时间周期。

图6 时间周期示意图Fig. 6 Time period diagram

1.2.3 模型求解及边界条件

将三维地热井简化成了一维线单元，考虑井筒内流体沿井筒轴向的渗流传热，井筒内流体与周围岩体的换热过程采用等效换热系数近似考虑，该系数考虑了井径、套管和水泥砂浆的导热性、流体性质和流量的影响。盖层传热过程可用热传导方程：

模型初始条件和边界条件：热储初始温度60 ℃，井口流速0.09 m/s，回灌温度20 ℃，储层初始水压为静水压力，井筒内初始水压也按静水压力考虑，运行开始后(t＞ 0)，井口设为流速边界条件，井筒内初始温度等于围岩初始温度。

模型的初始温度场可根据地温梯度 ΔT来确定，顶面温度边界条件设为全年平均气温，侧面温度开边界和水位边界，底边温度开边界和不透水边界，如图7 所示。模型所用热储岩石热物理参数如表1 所示。

图7 模型边界条件示意图Fig. 7 Model boundary conditions

2 模型验证及当前模式分析

2.1 计算模型与参数

根据不同地层的增温率，Surfer 软件插值得到不同地层的地温梯度分布，雾迷山组下伏岩层地温梯度数据未掌握，根据前人研究基础，取1.5 ℃/100 m。选择6 眼温度监测井进行拟合，调整模型的热物性参数使得温度监测井中的模拟温度与出水温度保持一致。结果表示温度与监测温度吻合结果较好，图8 为温度监测井和温度拟合数据，计算得到模型的初始温度场分布如图9 所示。

图8 温度拟合曲线Fig. 8 The temperature fitting curves

图9 初始温度场分布 /(℃)Fig. 9 The initial temperature distribution

利用区内5 眼地热井的抽水试验数据进行水位拟合，用裘布依公式和吉哈尔特降压影响半径经验公式迭代计算出监测井附近渗透系数的参考值。模型取大落程的渗透系数作为参考值，用Surfer软件插值得到储层的渗透系数分布，基于抽水试验(5 眼监测井)和回灌试验(京通-4)的抽水量、降深数据与水位监测数据进一步调整热储的渗透系数分布(表2)，最终使模拟水位与实测水位较好吻合，如图10 所示。

表2 校正前后渗透系数对比表 /(m/d)Table 2 Comparison of permeability coefficient before and after calibration

图10 抽水、回灌试验水位拟合曲线Fig. 10 Fitting curves of water level of exploitation and recharge tests

2.2 真实采灌模式评价

区域内共25 眼地热井，3 眼回灌井(通热灌-9、京通灌-4、京通灌-6)。结合地热井现有开采能力以及北京冬季供暖需求，假设回灌水温25 ℃，当单井采灌量为1000 m3/d、1500 m3/d 和2000 m3/d 时，计算开采100 年后热储的温度场和水位场的演化。

地热回灌可能导致热突破，将开采温度下降2 ℃定义为热突破，关注回灌井分布密集区域温度场的变化。距离通热灌-9 最近的井是通热-15，井间距为182 m，冷锋面很快到达了通热-15，当采灌量分别为1000 m3/d、1500 m3/d、2000 m3/d，对应的热突破时间为8 a、5 a、4 a，开采井中温度逐渐下降并稳定在38 ℃左右，如图11 所示。距离京通灌-6 最近的井是京通-3，在采灌量较大的情况下(2000 m3/d)出现了热突破，热突破时间为82 a。京通灌-4 位置相对孤立，距离最近的京通-5 的井间距超过2000 m，不受到回灌影响。

图11 开采温度变化曲线Fig. 11 The temperature variation curves

回灌井分布较稀疏的区域，回灌不足可能会带来水位场的变化。在长期开采的情况下，开采井集中区，形成大面积的抽水漏斗，并随着采灌量的增加而增大。

即使回灌井稀少，但对维持热储压力、稳定开采水位起到了相当大的作用。随着回灌量的增大、采灌井间距的减小，其作用也越明显。但是不合理的回灌井布置，可能会导致开采井过早发生热突破(如通热灌-9 导致通热-15 过早发生了热突破)，因此，合理的采灌井布置才是实现研究区地热资源长期开发利用的有效途径。

3 优化开采模型研究

针对目前采灌方案的不足，在100%回灌的前提下，结合工作区内的开采井、回灌井提出三种布井方案， “一采一灌”、“一采两灌”和“两采三灌”。

结合区域内地热井的开采能力及蓟县系雾迷山组热储开采经验，单井的采灌量分别考虑1000 m3/d、1500 m3/d 和2000 m3/d，采灌井间距分别考虑300 m、500 m 和800 m 三种工况，回灌温度为供暖尾水温度25 ℃。现有采灌模式下水位场分布如图12 所示。

图12 现有采灌模式下水位场分布 /mFig. 12 Water level distribution with current exploitation and recharging mode

“一采一灌”布井方案中，采灌量大于1500 m3/d时，冷锋面运移到了开采井处，开采井在78 a 发生了热突破，采灌量增大为2000 m3/d，热突破时间缩短至56 a；采灌井间距为500 m 和800 m 时，无论采灌量为1000 m3/d、1500 m3/d、2000 m3/d，回灌冷锋面均未运移到开采井处(图13)，开采井温度也没有下降，当采灌井间距大于500 m 时，开采井的热突破风险小。“一采两灌”布井方案中，由于回灌井布置在开采井下游位置，在三种采灌量下，回灌冷锋面均未运移到开采井处，未造成开采井发生热突破，开采井的热突破风险最小(图13)。“两采三灌”的布井方案中，位于下游的开采井受到回灌冷锋面的影响较大，热突破风险高。

图13 热储温度场分布 /(℃)Fig. 13 Temperature distribution of heat storage

在三种采灌量下，下游的开采井均发生了热突破，热突破时间分别为96 a、72 a 和58 a；位于上游的开采井在采灌量小于1500 m3/d 没有发生热突破，当采灌量增大到2000 m3/d 时在85 a 发生了热突破，位于下游的开采井热突破风险极高，位于上游的开采井也面临着一定程度的热突破风险。

从水位计算结果来看：不同布井方案下，热储的水位场分布有所不同，在开采井附近会形成低压区，回灌井附近形成高压区(图14)。在三种布井方案下，随着采灌量和井间距的增大，开采井的水位逐渐下降，回灌井的水位逐渐上升。采灌量的大小直接影响了采灌井水位的高低，而采灌井间距越小，回灌对稳定开采水位的作用越明显。

图14 热储水位场分布 /mFig. 14 Water level distribution of heat storage

在“一采一灌”布井方案中，当采灌量增大到2000 m3/d 时，回灌井的水位会超过0 m， “一采一灌”布井方案无法实现100%自然回灌。而在“一采两灌”和“两采三灌”布井方案中，由于回灌井数量的增加而减小了单眼回灌井的回灌量，使得回灌井水位均在0 m 以下，可实现100%自然回灌。

优化开采模式研究表明区内回灌井布设需要满足以下条件：

① 采灌井间距大于500 m，采灌量不超过2000 m3/d；

② 大地流场方向由西北流向东南，回灌井需布置在开采井下游；

③ “一采两灌”布井方案，部分开采井集中的地区可采用“两采三灌”；

④ 回灌井的布设位置不应对周围其他地热井造成干扰。

根据以上原则区内已有的距离其他井较远的开采井匹配回灌井时，采用“一采两灌”的布设方案(若不新增开采井，无法实现“两采三灌”)；部分距离较近的开采井，如京通-1 和通热-19、通热-9 和通热-15、京通-3 和京通6、园-1 和园-2 等采用“两采三灌”的布设方案，以节约运维成本。具体的回灌井布设方案如图15 所示，区内共布设回灌井46 眼。在该布井方案下，集中开采区的抽水漏斗范围得到了有效的控制，稳定了开采水位，在100 年的开采年限里，开采温度没有明显的下降，说明开采井未发生热突破(图16)。对比目前采灌模式下监测井的水位可知，优化采灌模式下能有效恢复并稳定开采水位不超过150 m(图17)，说明优化采灌模式能满足研究区内地热资源的长期开发利用。

图15 回灌井布设方案Fig. 15 The distribution plan of recharge wells

图16 优化采灌模式模拟结果示意图Fig. 16 The simulation results of optimized exploitation and recharging mode

图17 两种模式下监测井水位对比曲线Fig. 17 The curves of water levels monitored under the two modes

“一采两灌”布井方案，回灌井布置在下游可大幅降低开采井的热突破风险，且能满足100%自然回灌的需求，通过对采灌井间距的敏感性分析可知，最适合研究区的采灌井间距应大于500 m，通过对采灌量的敏感性分析可知，为确保开采水位不超过150 m，研究区的采灌量不宜超过2000 m3/d。优化采灌模式下能有效恢复并稳定开采水位不超过150 m，优化采灌模式能满足研究区内地热资源的长期开发利用。

4 结论

本文基于城市尺度地热群井高效数值模拟方法，建立了工作区三维热储数值模型，并通过对温度和水位监测数据的拟合，校正了模型参数，为后面地热资源优化开采模式研究奠定了基础。主要结论如下：

(1)通过对目前采灌模式的模拟研究，回灌对维持热储压力和稳定开采水位十分重要，并随着回灌量的增大以及采灌井间距的减小，作用越明显，不合理的回灌井布置，可能会导致开采井过早发生热突破，因此，合理的采灌井布置才是实现工作区地热资源长期开发利用的有效途径。

(2)通过研究不同布井方案下的井组模型得到了最适合工作区的布井方案为“一采两灌”和“两采三灌”，且当采灌井间距大于500 m，采灌量不超过2000 m3/d 时，能保证100%回灌且在100 年开采年限里，开采井不发生热突破，开采水位不超过150 m 的需求。提出了工作区内的优化采灌布井方案原则，在满足布设原则的优化采灌模式下，开采井不会发生热突破，且能有效恢复并稳定开采水位不超过150 m，实现了区内地热资源的长期开发利用。