谭志容，高宗军，赵季初

(1.山东省鲁北地质工程勘察院，山东 德州 253072；2.山东科技大学地科学院，山东 青岛 266590)



基于Modflow软件的德州市城区馆陶组热储数值模拟研究

谭志容1，高宗军2，赵季初1

(1.山东省鲁北地质工程勘察院，山东 德州 253072；2.山东科技大学地科学院，山东 青岛 266590)

德州市位于山东省西北部，区内地热资源丰富，城区已有地热井75眼，地热资源大量开发，对地热资源的管理提出了新的要求，本文以德州市城区为研究区，充分收集资料，采用Modflow软件建立德州市城区馆陶组热储模型，对馆陶组热储的埋藏分布规律及城区地热田的地热水储量进行计算评价，将馆陶组热储概化为平面上无限延伸且顶、底板为隔水层的均质含水层。采用非稳定流抽水试验求取的水文地质参数对模型进行识别，利用两眼长期监测孔的水位动态观测数据及历年来地热资源的开采情况对模型进行校正。校正过程中发现与热储砂岩互层分布的泥岩压密释水是热储地热水的主要来源之一。以地热水水位降幅2 m/a为限制条件，德城区馆陶组热储地热水的允许开采量为532万 m3/a。

德州市城区；馆陶组热储；数值模拟；通用水头边界；允许开采量

德州市位于山东省西北部，在大地构造单元上处于华北平原的东南部。自1996年在德州市城区开展地热资源普查、施工了两眼探采结合孔以来，区内地热资源开发蓬勃发展，并带动了周边地区地热资源的开发。目前德州市各县级城区均有数眼至数十眼地热井不等，德州市城区已有地热井72眼。

地热资源的大量开发，对地热资源的管理工作提出了新的要求，其中地热水可采资源量的确定是整个地热资源管理工作的核心。根据规范要求：“在地热田的勘查程度比较高，并且具有一定时间的开采历史，具有比较齐全的监测资料时，应建立地热田的数值模拟模型，用以计算/评价地热储量，并作为地热田管理的工具”[1]。

本文在充分收集整理以往勘探成果的基础上，借鉴前人在地下水数值模拟方面取得的成功经验[2-6]，采用Modflow软件建立德州市城区馆陶组热储的数值模型，对德州市城区地热田的地热水储量进行计算与评价。

1 地热地质条件

1.1 地质构造

德州市城区在大地构造单元上位于临清坳陷的东北部，东邻埕宁隆起、西接沧县隆起，北与黄骅坳陷相通，为两隆夹一坳的地堑型构造(图1)。其内沉积了巨厚的新生代地层，是地下热水的赋存场所，其中馆陶组下段砂岩是区内地热资源开发的目标热储[7-8]。

图1 模拟区构造位置纲要图

1.2 馆陶组热储特征

馆陶组热储在整个华北平原内皆有分布[9]，具有下粗上细的正旋回沉积特征，底部砂砾岩层分布较稳定。平面上热储呈层状连续分布，垂向上热储与泥岩呈互层分布。热储岩性主要为河流相、冲积扇相的细砂岩、粗砂岩、含砾砂岩与砂砾岩，砾石呈半圆状，磨圆度中等，矿物成分主要为石英与长石。

模拟区馆陶组热储的厚度与埋深受基底构造的影响，底板埋深1 300～1 600 m，热储厚120～200 m(图2)，单井涌水量80～100 m3/h，水化学类型为Cl—Na型，矿化度5 g/L左右。

图2 模拟区馆陶组热储厚度与底板埋深图

2 基于Modflow软件的数值模型

2.1 概念模型

根据区内馆陶组热储特征，为降低模型的复杂程度，本次将整个研究区概化为一个均质各向同性的稳定含水层。热储顶底板均分布有厚度巨大的泥岩层，在模拟期内可视为隔水层。热储在整个华北平原内均有分布，相对于整个华北平原而言，可将热储概化为平面上无限延伸的含水层，其概念模型如图3所示。

图3 热储概念模型图

2.2 无限含水层边界的概化

2.2.1 边界条件的分类

地下水流动问题的边界条件主要分为给定水头边界条件(第一类边界)、给定流量边界条件(第二类边界)与潜水面边界(陈崇希等.1999)[10]。

Modeflow 软件将模拟区划分为活动单元(Active Cells)与非活动单元(Inactive Cells)，非活动单元与活动单元之间的边界可认为是零流量边界(特殊的第二类边界)。此外，该软件还给出了8种不同的含水层边界条件，其中适宜本概念模型边界条件的有两种，即定水头边界(Constant Head)与通用水头边界(General Head)。

定水头边界(第一类边界)可通过边界附近的观测井监测资料，给出不同时间段的初始时水位与结束时水位，软件自动插值确定各模拟期边界的水位，当模拟区边界附近有较多的监测井时，这类边界条件能较真实地反映模拟区与外界的水量交换。

通用水头边界(第二类边界)通过计算该边界处活动单元的水位与模拟区外一定距离处给定的水位之间的差值，计算流入或流出模拟区的水量。

2.2.2 模型边界的选择

根据概念模型，热储概化为四周向外无限延伸的均质各向同性含水层，最为适宜的边界条件应该采用给定水头边界，但由于模拟区边界附近无监测井分布，没有办法给定各模拟时间段热储的水位。若采用零流量边界来定义模拟区四周的边界，则模型的范围要外扩，至少达到模拟区抽水井的影响范围之外，采用这种方法不能重点突出模拟区的水位变动情况。

因此，本次采用通用水头边界来刻画这种无限延伸的含水层边界条件，即假定离模拟区边界一定距离处的水位保持不变，这种假定也与实际情况基本相符。

2.3 模型的建立与识别

2.3.1 单元剖分

本次以德州市城区的地理底图作为模型建立的底图，将模拟区概化为一个南北长37 453 m，东西宽35758m的矩形区域。

根据地热井的钻探与成井资料，垂向上将整个模拟区概化为深度0～1 300 m、1 300～1 550 m及1 550～2 000 m三个地质层。其中1 300～1 550 m为区内地热井的取水层，即热储层，属馆陶组下段，地层岩性为砂砾岩、粗砂岩及中粗砂岩与泥岩的互层，其中砂岩占地层厚度的比例60%左右。0～1 300 m为热储盖层，从上至下为第四系、明化镇组及馆陶组上段；1 550～2 000 m为热储下伏层，从上至下为东营组与沙河街组。

平面上将整个模拟区划分为100行×100列的均匀网。

2.3.2 抽水井数据处理

德州市城区地热资源开采始于1998年，截止目前共有地热井72眼，其中1998年投入使用2眼，2002年新增3眼；2003年新增7眼，2004年新增4眼；2005年新增11眼；2006年新增24眼；2007年新增3眼；2008年新增8眼；2009年新增6眼；2010年新增1眼；2011新增3眼。

各地热井的集中开采时间段为11月15日～至翌年3月15日，约120 d，单井开采量约80 m3/h；其余时间段开采量极小，在模型中可以忽略。

因此，各井自投入使用开始，每年开采应力期概化为两段：即每年的11月15日～至翌年3月15日，单井开采量1 920 m3/d；每年的3月16日～11月14日开采量为零。以1998年1月1日为模拟的初始时间，以2016年2月21日为模拟的结束时间，总时长为6 625 d。

2.3.3 水文地质参数及边界

1) 渗透系数(K)与贮水率(Ss)

根据非稳定流抽水试验结果，区内馆陶热储的导水系数为944.5 m2/d，贮水系数(弹性释水系数)为2.48×10-4[11]。本次模型的热储厚度取值250 m，根据导水系数公式：

T=K·M

式中：T为导水系数(m2/d)；K为渗透系数(m/d)；M为热储厚度(m)。

经计算可得，热储的平均渗透系数为3.778 m/d。

根据贮水率(Ss)的定义，当含水层水头变化一个单位时，从单位体积含水层中，由于水体积膨胀(或压缩)以及介质骨架的压缩(或伸长)而释放(或贮存)的弹性水量[12]。根据贮水系数公式：

μe=Ss·M

式中：μe为贮水系数(无量纲)；Ss为贮水率(1/m)。

经计算，热储的平均贮水率为9.992×10-7(1/m)。

2)边界的输入

本次假定在抽水影响范围以外的热储保持定水头不变，其关键是影响范围的确定，采用裘布依公式进行估算[13]：

式中ΔS为距开采中心R处与开采中心的热储水位差(m)；Q为地热水开采量，本次取72眼井的共同开采量(m3/d)；R为影响半径(m)；r为72眼井圈定的假想圆半径(m)；

根据模拟区两观测孔的最大水位降约68 m左右，推算的R值为135 km左右，即将通用水头在离模拟区外135 km处的水位设为初始值“零”。

2.4 模型的校正

利用模拟区两眼观测孔长期水位动态监测数据对模型进行校正。两眼监测孔的编号分别为DR14-1与DR34-1。运行模型后，结果显示，计算值与观测值的拟合度差，模型参数不合理，需要校正。

通过对各水文地质参数进行分析可知，由于抽水试验取得的弹性系数主要反应含水层的贮水/释水能力，而实际情况是热储为砂、泥岩互层，在长期抽水作用下，泥岩的压密释水不容忽视，因此贮水率值要增大；渗透系数主要是反映热储的渗透能力，区内热储较厚，泥岩对热储渗透性能的影响程度略轻，因此热储平均渗透系数也要适当增大；另外，受复杂的地质条件影响，侧向迳流补给也没有想象中那么大，在通用水头边界赋值时，在不改变距离R值时应将模拟区外热储的渗透系数减小(相当于增加了R的距离)。

根据以上分析，经反复调试，最终确定校正前后的水文地质参数如下表所示(表1)：

表1 校正前后水文地质参数对比表

注：K’为通用水头边界赋值的模拟区外热储渗透系数

经参数校正后，检查校正结果图(Calibration graph)发现，计算值=观测值的直线位于95%置信区间的中心，绝大部分观测值位于95%区间内(图4)，说明模型已校正。

另外，两个观测井的水位降深观测值与计算值的拟合度较高(图5)，符合《地下水资源管理模型工作要求》(GB/T 14497-93)中关于“水位拟合时相对误差小于10%的节点必须占已知水位节点的70%以上”[14]的规定。

图4 模型校正前计算值与实测值校正图

图5 DR14-1、DR34-1两眼观测孔水位计算值与观测值拟合曲线

3 地热水允许开采量计算结果

经校正后，所建立的数值计算模型具有较好的准确度。利用该模型既可用来预测不同开采条件下的水位变动情况，亦可通过控制热储水位的年降幅，求得该允许的水位年降幅条件下的开采量作为地热水的可采资源量。参照规范[1]，以年压力下降速度不大于0.02 Mpa(相当于2 m/a)作为限制条件，运行该模型，求得德州市城区内馆陶组热储地热水的可采资源量为532万 m3/a。

4 结语

(1)德州市城区馆陶组热储地热水开发始于1998年，目前已有地热水井72眼。地热水机制开采期是每年11月15日至来年的3月15日，单井开采量为80 m3/h。

(2)由区域地质条件，德州市城区馆陶组热储可概化为厚度250 m、顶底板为隔水层的平面无限延伸含水层，模型采用了通用水头边界。

(3)分析研究发现，开采条件下热储的贮水率大于单井抽水试验时求取的贮水率，说明泥岩压密释水是热储地热水的重要来源，故贮水率上调了1个数量级。由此可见，长期地热水开采，势必造成不可恢复的地面沉降，应引起重视。

(4)以地热水位降幅2 m/a为限制条件，德州市城区馆陶组热储地热水可采资源量为532万 m3/a。

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

[2]薛禹群.中国地下水数值模拟的现状与展望[J].高校地质学报.2010,16(1):1-6.

[3]孙从军,韩振波,赵振，等.地下水数值模拟的研究与应用进展[J].环境工程.2013,31(5):9-13.

[4]武强，朱斌，徐华，等.MODFLOW 在淮北地下水数值模拟中的应用[J].辽宁工程技术大学学报.2005(4):503.

[5]沈媛媛,蒋云钟,雷晓辉，等.地下水数值模拟中人为边界的处理方法研究[J].水文地质工程地质.2008(6).

[6]童彦钊, 吴平, 韩强强. 地下水数值模拟软件简述[J]. 科技信息.2012, (29):48-48.

[7]颜世强,潘懋,邹祖光，等.山东德州凹陷地下热水地球化学特征及成因[J].中国地质.2007,34(1):149-152.

[8]赵季初.山东省德州市城区馆陶组热储地热资源评价[D].北京:中国地质大学(北京).2007.

[9]陈墨香.华北地热[M].北京：科学出版社.1988，90-91.

[10]陈崇希,林敏.地下水动力学[M].武汉:中国地质大学出版社(武汉市洪山区鲁磨路31号).1999:32.

[11]赵季初,谭志容,王成明.德城区馆陶组热储水文地质参数探讨[J].地热能.2008(2):16-18.

[12]张惠昌. 关于贮水系数与贮水率意义的商榷[J]. 勘察科学技术.1987(4).

[13]JGJ/T111-98，建筑与市政降水工程技术规范[S].23-24.

[14]GB/T 14497-1993，地下水资源管理模型工作要求[S].

Study on Numerical Simulation of thermal reservoir of Guantao formation in Dezhou based on Modflow software

TAN Zhi-rong1，GAO Zong-Jun2，Zhao Ji-chu1

(1.Shandong Lubei Geological Engineering Investigation Institute, Dezhou, Shandong, 253072, China；2.Academy of Sciences of Shandong University of Science and Technology, Qingdao, Shandong, 266590, China)

Shandong province is located in the northwest region, rich in geothermal resources, the city has 75 geothermal wells, geothermal resources development, put forward new requirements for the management of geothermal resources, this paper takes Dezhou city as the study area, to collect sufficient data, using Modflow software to establish the Dezhou City geothermal reservoir of Guantao formation model. Water distribution and urban buried geothermal field of reservoir of Guantao formation in the evaluation calculation, the reservoir of Guantao formation is generalized to plane and infinite top floor aquifuge homogeneous aquifer. The unsteady flow pumping test of hydro geological parameters obtained from the test is used to identify the model, with both eyes the long-term monitoring of water level dynamic observation data of holes and over the years, exploitation of geothermal resources of the model. The sandstone and reservoir interblended mudstone compaction released water is the main reservoir of geothermal water found in the process of correction One of the sources is to reduce the geothermal water level of 2m/a as the limiting condition, the allowable exploitation of geothermal water in the Guantao formation of area is 5 million 320 thousand

Dezhou city; guantao reservoir; numerical simulation; general head boundary; allowable exploitation

2017-02-20

谭志容(1975-)，女，湖北荆州人，高级工程师，主要从事水文地质方面工作。

高宗军(1964-)，男，山东泰安人，教授，研究方向：水文地质。

P314.1

A

1004-1184(2017)03-0019-03