天津市地热资源可持续利用管理模型研究

刘杰，宋美钰，于彦，边宗斌，康楠，秦莉红

（1.天津地热勘查开发设计院，天津300250；2.天津地质资料馆，天津30000）

天津市地热资源广泛应用于供暖、洗浴、生活和工农业用水等诸多领域。自上世纪八十年代以来，地热资源的开发利用进入了快速发展时期，开发利用形式也从粗放型逐渐发展到集约化阶段。面对地热资源需求的快速增长和开采规模的不断扩大，如何维持可持续开发利用，为地热资源最优化开发利用的采灌布局提供依据，是目前我市地热行业急需解决的问题。本文以天津潘庄凸起区为例，说明应用优化管理模型进行地热资源优化管理。

1 模型建立及识别

1.1 水文地质模型

选取Ⅳ级构造单元潘庄凸起区为本次的研究区，具体边界：北部至汉沽断裂，西部以天津断裂为界，东部边界为沧东断裂[1]，南部为自然流量边界，总面积374 km2（图1）。研究区内主要以基岩裂隙型蓟县系雾迷山组热储层为主要开采层位，共有地热井18眼。

研究区内雾迷山组热储层普遍分布，岩性以灰褐色含灰泥晶白云岩、泥晶白云岩为主，含石英粉砂泥晶砂屑。钻井揭露顶板埋深为1 752～2 016 m，揭露厚度为480～1 032 m，揭露的最大深度3 634 m，尚无钻孔钻穿，故为了模型计算，热储层底板统一以4 000 m计算。

根据本区的地热地质、水文地质条件和地热流体开发利用特点、同位素分析结果：计算区西侧天津断裂为隔水边界、但不隔热；北部汉沽断裂为定流量流入边界；东部沧东断裂为热流体导水导热通道；南部为自然流量边界[2]。

1.2 数学模型

对所研究的流体系统封闭边界Гn包围的任意子域Vn积分。TOUGH2求解的质量和能量平衡方程的一般形式可以写为：

式中：MK代表单位体积的质量或能量（κ=1、2、…… NK，标记物质组分水、空气、CO2、溶质、……）；φ代表孔隙度；sβ代表状态β的饱和度（即β态占据的毛孔体积部分）；ρβ代表状态β的密度（kg/代表组分κ出现在状态β中质量比；F代表表示质量或能量通量（流体流量）；q代表源汇项；n代表表面单元dГn上的指向dVn的法向矢量。同时根据Darcy定律多相形式给出各状态的通量：

式中：ρβ代表状态β的密度(kg/m3)；υβ代表状态β的Darcy速率（体积通量）（m/s）；k代表绝对渗透率（mD）；krβ代表状态β的相对渗透率（mD）；μβ代表粘度（kg/m/s）；Pβ代表状态β中的流体压力，它是相对其他状态（通常参照气态）压力P和毛细压力Pcβ（≤0）的总和；g代表重力加速度向量。

1.3 数学模型识别

1.3.1 研究区剖分、初始条件及时间离散

根据雾迷山组热储层的埋藏分布条件、地质构造分布特征、富水性特征及地热井的布局和水位动态监测情况[3]，利用TOUGH2.0模拟软件将研究区剖分为4 580个单元（图2）。

研究区有较系统的水位常观资料始于2002年，所以以2002年1月为模型模拟的初始时刻。

拟合时段分供暖期（120天）与非供暖期（245天），时间步长为30天，模型识别时间到2014年10月。

1.3.2 边界补给量

边界的具体补给条件根据水位等值线图及历史水位拟合确定，相应的单宽补给量可根据水位等值线图计算，并在拟合过程中调整（表1）。

图2 研究区剖分图Fig.2 Subdivision graph maps in the study area

1.3.3 参数介质分区

计算水文地质参数选用是依据三维模型的计算需要，主要采取分区赋值的方法。根据雾迷山组热储层的埋藏分布和控热构造分布特征，利用所收集的降压试验资料进行综合统计分析，按所求水文地质参数进行初始分区。此次共收集9眼雾迷山组地热井降压试验数据（表2），根据资料[4]将该组热储层水文地质参数划分为3个初始分区（图2）。

1.3.4 模型识别与验证

模型拟合的开始时间选择在2009年10月，模拟终止于2014年10月，划分为61个时间步长（每个月为一个时间步长）。

对雾迷山热储层的观测水位按液面20℃温度进行统一换算，利用模型计算出的水位值与观测值进行拟合对比（图3）。可以看出：计算值与观测值在变化趋势和数值都很接近，大部分实际观测值点都在模拟曲线上[5]，计算值与观测值平均相对误差DL-19为1.38%；DL-19B为2.39%；DL-40为2.12%；DL-40B为1.98%，计算值和观测值较为一致。

表1 边界单宽补给量(m3/d·m)Tab.1 Single wide recharge boundary(m3/d·m)

表2 热储模型主要参数表Tab.2 Main parameter of the thermal reservoir model

从上述综合分析可以看出，研究区数值模拟模型所对应的地热地质模型与实际地热地质模型是基本吻合的，同时也说明该模型能基本上合理地反映示范区的地热地质条件，为下一步建立地热资源管理模型奠定基础。

2 管理模型的建立

研究区地热资源优化管理的目标[6]：一是在尽可能满足地热生产需求的条件下，实现全区水位降幅最小；二是针对目前的地热开发利用现状，尽量减少地热井之间相互干扰，充分发挥模型调控功能。

2.1 目标函数的建立

研究区地热资源管理面对的是整个地热资源系统，目的是规划、调整地热资源系统的时空人为输出，使地热流体输出趋近于既定目标[7]，同时使地热资源系统的天然补给朝着有利于实现系统功能的方向发展以使有限的地热资源发挥最大的经济、环境和社会效益。

故取示范区中地热井的开采/回灌量为模型的决策变量，即：该模型的目标函数为使n年的规划阶段、m个地热井的地热资源总开采量最大：

式中：maxZ代表最大地热资源总开采量；m代表地热井数量；n代表规划阶段（年）；Q代表地热井开采量。

本次规划阶段分为30年，其中计算A1阶段（1年）、A5阶段（5年）和A30阶段（30年），共有12眼地热井参与计算，即可得总开采量为：

式中：maxZ代表最大地热资源总开采量；Q代表地热井开采量。

2.2 约束条件的建立

2.2.1 约束条件的选择

约束条件首要的是水位降深约束[8]，要求各井的水位降深不能超过某个给定的值。如果是给定的值，该值要考虑三个因素：一是不应超过地热规划中规定的热储层中水位最大埋深；二是经济效益，开采井降深过大会导致能耗过大；三是不应超过抽水试验达到的最大降深。

图3 研究区地热井水位埋深（20℃）计算值与实测值拟合图Fig.3 Water depth(20℃)calculated values and measured values fitting figure of the geothermal wells in the study area

式中：β（i，j，t）代表单位脉冲响应函数，表示由在第j点脉冲的作用，在第t时段末第i点的水位降深响应；Q（i，j）代表第t时刻作用于j点施加的脉冲量(开采/回灌量)；Si代表第i点的水位埋深约束。

②水位降深极限约束：

式中：S(i，t)代表第t时刻i点的水位降深；S0(i)代表第i点的水位降深约束。

③决策变量上下限约束：

④非负性约束，表示各决策变量为非负：区内的地热资源开发利用更趋于采灌平衡状态（控制每年降幅不大于1 m）。

式中：Q(j，t)代表第t时刻作用于j点施加的脉冲量（开采/回灌量）；S(i，t)代表第t时刻i点的水位降深。

2.2.2 约束条件的具体数值

（1）水位埋深约束

根据天津市地热资源管理办法，雾迷山组热储水位埋深总体不超过200 m。

（2）决策变量上下限约束

决策变量上下限约束是指在满足生产需求的前提下，对开采量进行合理的上下限控制。针对示范区内12眼参与计算的地热井进行调查，根据近三年的开发利用情况制定约束数值（表3）。

（3）水位降深约束

水位降深约束条件指在分布于研究区内的18个水位降深约束点处，把每个约束点的水位降深约束在一定水平上的一组数据。以基本维持目前水位为基准，不超过控制开采区标准为要求，并且分别计算水位降深为0 m/a，1 m/a，3 m/a的优化方案。

为了下文方便描述，将水位降深约束为0 m/a定为优化方案一、水位降深约束为1 m/a定为优化方案二、水位降深约束为3 m/a定为优化方案三。为了体现优化管理模型在地热资源管理上的管理能力，另外选取一种方案（优化方案四），即保持示范区内开采总量和回灌总量不变的情况下，通过管理模型调整各个地热井的开采量或回灌量，进而让整个示范

表3 决策变量约束统计表Tab.3 Statistical data of the decision variables constraints

3 计算方法与程序设计

3.1 计算方法

本项目采用Matlab编程进行优化计算[9]，对管理模型进行最优求解。利用Matlab编程来进行线性优化，采用Matlab中单纯形法（Linprog函数）进行求解[10]。

3.2 程序设计

现以第一时段（即2015年管理模型）为例，以下为Matlab程序设计：

Linprog线性规划函数

x=Linprog（f，A，b）解如下形式的线性规划问题

minfx

x

Subject to：Ax≤b

x=Linprog（f，A，b，vlb，vub）参数vlb、vub给出设计变量（开采/回灌量）的上下边界约束，即vlb≤x≤vub。

x=Linprog（f，A，b，vlb，vub，x0）设置初始值

[x,fval,exitflag,output,lambda]=linprog（c,a,b,[],[],vlb,vub）

返回目标函数最优解x，和在x处的值：fval=fT*x，是否存在exitflag标志，优化解结构体output，拉格朗日乘子结构体lambda。

4 优化结果及模型检验

4.1 优化结果

对12眼地热井不同条件下各时段分别建立MATLAB公式进行求解，结果见表4。

表4 不同优化方式的总量（×104m3）Tab.4 The quantity of different ways of the optimization（×104m3）

从表4可以看出：1）2014年蓟县系雾迷山组热储层年开采量为213.22×104m3，各时段的开采量基本能满足示范区地热资源的开发利用及管理模型建立的限制条件。2）在该区域每年降深为0 m或1 m时，优化管理后回灌量都比开采量大，只有当降深为3 m的时候，优化管理后的开采量比回灌量大，证明该区域地热资源开发利用相对密集，地热井相互之间影响较大，必须增加回灌来保证水位降深不断变小。3）从表中可以看出水位降深越大，优化开采量越大；在同一开采方式下，年份越长，优化开采量也越大。所以可以利用管理模型对区域内地热资源进行资源调控，对降深大的区域按照优化管理后的开采量进行控制，把优化开采量定为年开采指标。

4.2 模型检验

根据优化方案四（以2014年开采/回灌量为依据，1 m/a水位降深为约束）得出的优化方案以及现有地热资源开发利用情况，基于建立好的TOUGH2.0数值模拟平台对雾迷山组热储层压力场动态变化情况进行预测，具体预测数值如表5。

分别从优化方案四中取出相同的四眼地热井（DL-19、DL-19B、DL-40、DL-40B），与按2014年开采方案开采预测值进行对比（图4）。

从上述数据可以看出：按优化方案四（水位降幅约束为1 m，按2014年开采/回灌量）计算的30年优化方案代入数值模型中，出来的结果每眼地热井平均降深最大为1.8 m，最小为0.6 m。

图4 优化方案四与现用方案对比曲线图Fig.4 Comparison chart of the optimized scheme four and current solutions

综合分析，将三个优化方案结果分别代入数值模型进行模拟，得出每眼地热井实际地热流体水位降深与优化方案前提降深值基本吻合，说明该管理模型能合理计算出优化开采量，比维持现状开采更为合理。

表5 按优化方案四计算的各年份数值模拟结果Tab.5 The numerical simulation from the optimized scheme four calculating results in each year

5 结论

天津地热资源优化管理模型的研究表明：面对天津地热资源需求的快速增长和开采规模的不断扩大，对其进行优化管理是使其能够可持续利用的有效方法。

科学调度合理开采地热资源，天津地区地热资源丰富，但如何可持续开采面临很大困难。在开采地热流体保证环境安全的同时，严格控制地区水位埋深，使地热资源开发利用达到“采补平衡”。今后地热资源开发利用都按本研究对各区的约束条件进行控制，在对约束条件论证后管理部门应考虑制定有关条例规定，保证天津市地热资源更精细化管理。

[1]天津市地质矿产局.天津市区域地质志[M].北京：地质出版社，1992.

[2]GB/T 11615-2010地热资源地质勘查规范[S].

[3]胡燕，林黎，等.天津市地热资源可持续开发潜力评价报告[R].天津:天津地热勘查开发设计院，2015.

[4]宗振海，刘杰，等.2015年天津地热资源开发利用动态监测项目[R].天津:天津地热勘查开发设计院，2015.

[5]宁立波，董少刚，马传明，等.地下水数值模拟的理论与实践[M].武汉：中国地质大学出版社，2010.

[6]国家技术监督局.地下水资源管理模型工作要求(GB/T 14497-93)[Z].1993.

[7]舒艳，王红旗，等.地下水管理模型研究进展[J].水文地质工程地质，2005，32（6）:85-90.

[8]卢文喜.地下水系统的模拟预测和优化管理[M].北京：科学出版社，1999.

[9]王沫然.MATLAB与科学计算[M].北京：电子工业出版社，2005，第2版.

[10]尹尚先.地下水管理模型中响应函数解法探讨[J].水文地质工程地质，2001，29（2）：33-35.