MODFLOW在天津滨海新区地热流体数值模拟中的应用

张芬娜，贾志，2，高亮，2，刘九龙

（1.天津地热勘查开发设计院，天津 300250；2.中国地质大学（北京），北京 100083）

MODFLOW在天津滨海新区地热流体数值模拟中的应用

张芬娜1，贾志1，2，高亮1，2，刘九龙1

（1.天津地热勘查开发设计院，天津 300250；2.中国地质大学（北京），北京 100083）

为了合理开发利用地热资源，本文以天津滨海新区桥沽-看财庄新近系明化镇组地热流体为例，应用Visual MODFLOW软件对明化镇组地热流体进行数值模拟。对模型的识别与检验表明，所建模型能够较好的反映明化镇组地热流体的地下水位动态特征。计算得出明化镇组100年水位接近150 m的地热流体可采资源量为10.7×106m3/a。通过预测按设计开采量进行开采后明化镇组地热流体地下水位的变化趋势，每天新增1000 m3开采量是可行的。

天津滨海新区；地热流体；数值模拟；MODFLOW

近几十年来，地下水模型已被广泛应用于地下水资源评价、预报和管理之中[1]。目前国内常用的模拟地下水位动态预测的软件有Visual MODFLOW、FEFLOW和GMS[2]等。有许多学者建立了MODFLOW或者FEFLOW地下水流模型，用水文观测井的地下水动态观测资料，对数值模型进行了识别和校验后，用于预测不同开采方案的地下水动态变化[3-8]。天津市地热流体地下水位模拟研究成果有胡燕建立的天津市滨海新区新近系馆陶组地热流体流场三维数值模型[9]和田光辉建立的灰色系统GM(1,1)模型[10]等。但在地热流体地下水位动态预测的模拟中利用Visual MODFLOW和FEFLOW软件的较少，Visual MODFLOW软件在计算模型的建立、计算结果的可视化、数据前后处理以及与其他软件数据信息交互等方面能力较好[11-14]。因此，本文运用Visual MODFLOW建立明化镇组地热流体三维水流模拟模型，计算明化镇组的可采资源量，并通过模拟地热流体地下水位动态，来判断增加设计开采量是否可行。

1 研究区概况

研究区地处天津滨海新区北部桥沽－看财庄一带，位于Ⅲ级构造单元黄骅坳陷的北部，横跨宁河凸起和北塘凹陷两个Ⅳ级构造单元[15-16]，区内及其附近主要发育的断裂有：汉沽断裂、宁河断裂、大神堂断裂、茶淀断裂和北塘东断裂（图1）。在当今经济技术条件下区内可开采的热储层有新近系明化镇组、馆陶组和古近系东营组[17]。本文主要以新近系明化镇组作为研究目的层。

2 热储特征

图1 研究区构造简图Fig.1 Sketch map of the structures

明化镇组作为一套曲流相碎屑岩，在全区均有

分布，储集条件良好，厚度发育稳定，具有北部宁河凸起浅，南部北塘凹陷深的特点[18]。顶板埋深376～506 m，砂岩占总厚度的25.7%～29%，单井涌水量54～85 m3/h，液面温度25～53℃，总矿化度439.1～584.3 mg/L，pH值8.15～8.69，水质类型以HCO3-Na型为主[19]，局部地区为HCO3·SO4-Na型水[20]。2012年明化镇组年开采量69×104m3，水位58～62 m，年降幅1 m左右[21]。

3 概念模型

计算范围：为了在地下水流数值模拟过程中更加准确地给定边界水位（或流量），在确定模拟区范围时，将模拟边界适当向外扩展到具有长观孔游乐港1（观3）的区域，将东南部边界扩展到海岸线，整个数值模拟区域总面积为423.4 km2（图1）。

边界条件设定：把明化镇组顶底板定为隔水边界。通过本区的地热地质、水文地质条件和地热流体开发利用特点、同位素分析结果可知，计算区新近系明化镇组热储层为层状结构，分布范围广阔，不是完整的水文地质单元。因此，侧向边界设为自然流量边界（通用水头边界）。

4 数学模型的建立

模拟区内水文地质参数随空间变化，参数在水平方向没有明显的方向性，因此地热流体模型可概化为三维、非均质水平各向同性的非稳定热流体系统，在笛卡尔坐标系下可描述为[22]：

式中：Ω为渗流区域；h为热流体的水位标高(m)；Kx、Ky、Kz分别为热流体储积层中x、y、z方向的渗透系数(m/d)；Kn为边界面法向方向热流体储积层渗透系数(m/d)；Kz为边界内侧热流体储积层垂向渗透系数(m/d)；μs为单位储水率（1/m）；σ为通用水头边界的阻力系数，σ=L/K；L为模型边界到通用水头边界的水平距离(m)；K为模型边界到通用水头边界之间热流体储积层平均渗透系数(m/d)；ε为源汇项(m/ d)；h0为初始水位分布(m)，h0=h0(x,y,z)；hb为通用水头边界上的水位分布(m)；n为边界面的法线方向；z→为热储层顶、底边界面的法线方向；Γ1为通用水头边界；Γ2为渗流区域顶、底边界。

5 模型数值模拟

5.1 时空离散

模拟区面积约514 km2。利用Visual Modflow程序对明化镇组水位动态进行模拟，采用矩形等距剖分，剖分成80×60个单元格，其中活动单元格3955个，非活动单元格845个。有效活动单元计算区面积423.4 km2，有效活动单元计算区即为本次工作的模拟计算区。

5.2 数据处理

以2012年4月明化镇组热储层动态长观资料为基础，对地热流体的温度、密度统一校正（热储层温度统一校正到40℃），利用Kriging优化插值法取得初始水位，通过模拟运算，修改流场，得到该层的初始流场。

源汇项包括侧向径流补给、层间越流补给、人工回灌量、人工开采等。将源汇项输入时，有具体井位的按坐标输入，其它单元的源汇项概化为井的形式，平均赋于区内的各个网格点中。其中层间越流量由Visual MODFLOW模型软件根据所赋相关参数自动进行处理。

5.3 模型的识别验证

本次模型识别验证采用试估－校正法，即反复调整参数，达到较为理想的拟合结果。本次模型识别时段选择2012年4月1日－2013年10月1日，明化镇组参与拟合的长观井共4眼（图1）。可进行地热流体水位历时曲线和区域流场的拟合（图2），计算出来的长观井地热流体水位动态变化与实际变化趋势一致，热流体流场拟合等值线与实测流场吻合较好（图3）。根据明化镇组热储层的埋藏分布和控热构造分布特征分为3区，各区的水文地质参数主要有：渗透系数K和弹性释水系数μ*。各区的水文地质参数初值主要参考

相关的报告①胡燕，林黎，孙宝成，等.《天津市滨海新区地热资源可持续开发潜力评价》天津：天津地热勘查开发设计院，2008.。运行Visual MODFLOW，可得到在给定水文地质参数和源汇项条件下的地下水位时空分布，通过拟合同时期热储层的流场和典型孔的动态过程曲线，识别热储层的参数，识别后的水文地质参数（表1）符合实际的水文地质条件。

图2 长观井水位埋深观测值与计算值曲线对比图Fig.2 Comparison between the water-level-burial-depth curve of the long-term-observation well and the calculated curve

图3 明化镇组热流体水位拟合图Fig.3 Minghuazhen geothermal fluid level fitting chart

表1 模型识别后水文地质参数取值Table 1 Values of hydrogeological parameter model identification

6 模拟模型的应用

6.1 地热可采资源量的计算

计算热流体可采资源量采用方法的不同，可导致结果差异较大。热流体可开采资源量多少取决于开采条件下的补给量，也与开采的技术经济条件和开采方案有关。在约束限制条件一定的情况下，开采井布局的合理性成为求取可采量大小的关键。虚拟井采用梅花型布局，集中布置在地热地质条件较好的地区，在地热地质条件较差的地区适当布置几眼虚拟井，代入模型进行计算。时间以2012年10月的地热流体流场作为开始时刻，2112年10月作为结束时间；开采量在计算中调节，求出100年地热流体水位最接近150 m的地热流体采量作为该热储层的可采资源量，结果见表2。

表2 明化镇组地热流体可采量表Table 2 Minghuazhen group geothermal fluid recoverable scale

6.2 开采预测

鉴于明化镇组热储层为模拟区主要开采层段，开采量较大，考虑到该热储层回灌效果不好的现状，本次预测暂不考虑回灌问题，故本次在现有开采井的基础上，在地热开采较少的地段增加少量开采虚拟井，本次布置4眼，开采量按1000 m3/d代入模型进行运算，以2012年10月为初始时刻，初始水位按60 m（20℃）计算，分别预测10年和30年地热流体水位的

变化情况（图4、图5）。从预测结果可以看出，30年后模拟区明化镇漏斗中心（汉沽城区）水位将达到96 m左右，年降幅约1.2 m/a。由此可见，在现有开采量的基础上，合理布设井位，每天新增1000 m3的开采量是可行的。

图4 明化镇组10年预测地热流体水位等值线图Fig.4 Minghuazhen group 10 year forecast geothermal fluid level contour map

图5 明化镇组30年预测地热流体水位等值线图Fig.5 Minghuazhen group 30 year forecast geothermal fluid level contour map

7 结论

（1）本文在水文地质条件和地热地质条件分析的基础上，运用Visual Modflow三维模型，计算出研究区100年明化镇组地热流体水位接近150 m的地热流体可采资源量为10.7×108m3/100 a。

（2）综合考虑数值模拟的预报结果以及区域以往动态资料认为：本区在原有开采量的基础上，合理布设开采井，每天新增1000 m3的开采量是可行的。

[1]Anderson M P,Woessner W W.Applied groundwater modeling simulation of flow and advective transport[M]. San Diego：Academic Press,1992,215-286.

[2]郭晓东,田辉,张梅桂，等.我国地下水数值模拟软件应用进展[J].地下水,2010，32(4)：5-7.

[3]李树文，王文静，崔新玲.基于Visual Modflow的成安县地下水动态研究[J].河北工程大学学报（自然科学版），2008，（4）：49-52.

[4]任印国，柳华武，李明良，等.石家庄市东部平原FEFLOW地下水数值模拟与研究[J].水文，2009，29（5）：59-62.

[5]贾金生，田冰，刘昌明.Visual MODFLOW在地下水模拟中的应用—以河北省栾城县为例.[J].河北农业大学学报，2003，26(2)：71-78.

[6]高卫东，孟磊，张海荣，等.FEFLOW软件在地下水动态预测中的应用[J].上海地质，2008，4：10-13.

[7]邵景力，崔亚莉，赵云章，等.黄河下游影响带（河南段）三维地下水流数值模拟模型及其应用[J].吉林大学学报（地球科学版），2003，33(1)：53-55.

[8]武强，朱斌，徐华，等.MODFLOW在淮北地下水数值模拟中的应用[J].辽宁工程技术大学学报，2005，24（4）：500-503.

[9]胡燕，沈键，刘杰，等.天津市滨海新区新近系馆陶组地热流体流场三维数值模拟方法[C].第十三届中国科协年会第14分会场-地热能开发利用与低碳经济研讨会论文集[M].北京:中国科学技术协会，2011：106-113.

[10]田光辉，曾梅香，程万庆，等.灰色系统GM（1，1）模型在天津地下热水水位预测中的应用[J].地下水，2009，31(6)：8-9.

[11]John Marler,Shemin Ge The Permeability of the Elkhom Fault Zone，South Park，Colorado[J]Groundwater, 2003，41(3):223-233.

[12]Wu Qiang,Yin Zhongmin,Wu Xun.Groundwater Study in the Weerselo Area,Overijssel,the Netherlands[M]. Beijing：Petroleum Industry Press,2001,47-74．

[13]Moench A F.Importance of the Vadose Zone in Analyses of Unconfined Aquifer Tests[J].Groundwater,2004,42 (2):223-233.

[14]张志忠，武强.河水与地下水耦合模型的建立与应用[J].辽宁工程技术大学学报，2004，23(4)：449-452.

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

[16]刘杰，宋美钰，田光辉.天津地热资源开发利用现状及可持续开发利用建议[J].地质调查与研究，2012，35(1):67-73.

[17]胡燕，高宝珠，靳宝珍，等.天津地热流体水化学分布特征及形成机理[J].地质调查与研究，2007，30(3):213-218.

[18]杨玉新，高亮，张芬娜，等.天津滨海新区桥沽-看财庄地热异常区普查报告[R].天津：天津地热勘查开发设计院，2009.

[19]潘家明，眭彪，李耀义.天津滨海新区中新生态城地热资源赋存条件探讨[J].地质调查与研究，2010，33(3):185-189.

[20]阮传侠，于彦，高宝珠，等.天津滨海新区地热流体水化学特征分析[J].地下水，2010，32(1)：51-53.

[21]田光辉，宋美钰，刘东林，等.2012年天津地热资源开发利用动态监测年报[R].天津：天津地热勘查开发设计院，2012.

[22]姬亚东，柴学周，刘其声，等.大区域地下水流数值模拟研究现状及存在问题[J].煤田地质与勘探，2009，37（5）：32-36.

Abstract:We carry out the deep-penetrating geochemical prospecting test for thick loess coverage area in the Yuncheng basin.This work adopts four methods simultaneously,including clay absorption phase metal ions(Cu Co Ni Zn etc.),thermal released carbonates(ΔC),conductivity(η%)and fluoride ion selective electrode.The test objects are two copper-nickel-polymetallic sulfide deposits buried in 200～300 meters deep under the loess.In the experiment,the four methods all captured the geochemical information from 200-meter-deep metal deposits covered by the thick loess in different degrees.It proves that the methods are effective and feasible preliminarily. There is obvious positive correlation between the anomaly intensity of the clay absorption phase metal ions(Cu Co Ni Zn)and the concentrating strength in the mineral ore.The indirect prospecting information including thermal released carbonates(ΔC)and conductivity(η%)shows the same tendency.Furthermore,we superimpose the anomaly information of the four methods and form the assemblage anomaly,which overcome the defect of single information weak and unstable,indicate the buried polymetallic ore effectively,and reflect the mineral ore anomaly more obviously.The cost of production and test of the four methods can be accepted by grass-roots unit,because the field sampling techniques are simple and easy in operation,consistent with conventional soil geochemical exploration methods.The development of the deep-penetrating geochemical prospecting method has the popularization value in similar soil coverage area.

Key words:thick loess;deep-penetrating geochemical method;clay absorption phase metal ions;thermal released carbonates(ΔC);conductivity;fluoride ion selective electrode;assemblage information

MODFLOW Application on the Numerical Simulation of Geothermal Fluid in Tianjin Binhai New Area

ZHANG Fen-na1,JIAZhi1,2,GAO Liang1,LIU Jiu-long1
（1.Tianjin geothermal exploration and designing institute,Tianjin,300250,China; 2.China University of Geosciences,Beijing,100083,China）

In order to rational development and utilization of geothermal resources in Tianjin Binhai NewArea, we take Minghuazhen geothermal fluid in the Qiaogu-Kancaizhuang Neogene Minghuazhen Formation as anexample,do numerical simulation with Visual MODFLOW software.The result Indicate that the calculation model can reflect the dynamic characteristics of groundwater level in the geothermal fluid.It is suggested that in 100 years,geothermal fluid geothermal fluid level close to 150 m recoverable reserves is about 10.7×106 m3/a. Accoding to the prediction of changes in the design of mining Minghuazhen geothermal fluid of underground water level,to add 1000 m3exploitation of dayis feasible.

Tianjin Binhai NewArea；geothermal fluid；numerical simulation；MODFLOW;calculation;Minghuazhen Formation

Study on Prospecting Test of Deep-Penetrating Geochemical Multi-Mathod in Thick Loess Coverage Area of the Zhongtiao Mountain

JIN Zhi-bin,ZHOU Xin-peng,ZHANG Shuang-kui,ZHANG Lu
(Geophysical and Geochemical Exploration Institute of Shanxi Province,Yuncheng,Shanxi,044000,China)

P314.1

A

1672－4135（2014）03－0237-05

2014-05-06

天津市国土资源与房屋管理局项目：天津滨海新区桥沽—看财庄地热异常区普查：（国土房任[2008]10号）

张芬娜（1981-），女，工程师，2005年毕业于中国地质大学（北京）水文与水资源工程，现主要从事地热方面的研究，E-mail:ying-ying22@163.com。