基于Visual MODFLOW的地下水数值模拟研究及水位预测

2021-03-18韩琳颜翠翠赵振伟吴婧马鹏飞

山东国土资源 2021年3期

韩琳，颜翠翠，赵振伟，吴婧，马鹏飞

(1.山东省地质矿产勘查开发局八〇一水文地质工程地质大队，山东济南 250014；2.山东省地勘局第二水文地质工程地质大队(山东省鲁北地质工程勘察院)，山东德州 253072)

0 引言

水是人类社会赖以生存和发展的不可代替的但可调控的重要自然资源，是构成区域生态环境系统的基本要素[1-2]。目前，济南市城市生活用水由地下水和地表水组成，占生活用水绝大部分的地表水的使用不仅要经过多道工序，而且水质、口感也与地下水有一定差距[3]。近年来，市政府下决心解决居民饮用优质地下水问题，已开发长孝水源地，设计供水能力10万m3/d，目前已实现供水4万m3/d。因此对长孝地区的水文地质条件开展研究及时且必要，但在地下水的开采使用过程中，往往存在对地下水分布与动态变化了解不清的盲目开采，因此了解地下水的分布与动态变化就显得格外重要。但地下水的储藏条件特殊，常规的探测方法不能及时全面地了解地下水的赋存状态[4-5]，而数值模拟方法因其直观全面的优点，非常适合用于地下水的研究中。地下水数值模拟的研究方面，不少研究者展开了深入研究。于翠翠等[6]应用GMS软件三维地下水流数值模拟模型，对济南泉域内岩溶地下水进行数值模拟和水平衡分析，并应用时间序列分析法对泉水水位动态进行了预测。朱君等[7]利用FEFLOW软件模拟计算了丘陵山区地下水流动特征下氚的迁移规律。不少研究者就地下水数值模拟方法进行了系统性总结，分析了各种软件的适用条件[8-15]。Visual MODFLOW软件在水位预测方面应用广泛，不少学者应用该软件对不同地点的水位进行了预测，预测结果达到了工程要求[16-20]。本文基于以上研究利用Visual MODFLOW软件对模型进行数值模拟计算，展开地下水水位预测相关研究。

1 研究区范围及水文地质背景

1.1 研究区边界条件

济南地区在大地构造上横跨新华夏第二隆起带的鲁西隆起及新华夏第二沉降带的鲁西北坳陷区。南部是北倾单斜构造，北部是松散堆积物的凹陷区。南部单斜构造中发育有多条规模较大的NNW向断裂，如牛角店断裂、马山断裂、千佛山断裂、文化桥断裂、东坞断裂、文祖断裂等,这些NNW向断裂自东向西大致等距分布，将单斜构造分割为若干个断块，这就形成了济南地区地质构造的基本特征。根据区内断裂构造的水理性质将济南地区岩溶地下水划分为5个水文地质单元(图1)。

1—第四系；2—石炭系；3—奥陶系；4—寒武系；5—泰山岩群；6—侵入岩；7—岩脉；8—奥陶系顶板埋深等值线(m)；9—水文地质单元分区界线；10—研究区范围界线图1 济南地区水文地质单元分区略图

长清-孝里铺水文地质单元是济南地区五大水文地质单元之一，面积935.40km2，各地质单元的面积见表1。市政府为解决居民饮用优质地下水问题开发的长孝水源地就位于长清-孝里铺水文地质单元内。长清-孝里铺水文地质单元边界条件及与周边水文地质单元水力联系如下：

表1 济南水文地质单元分区及分区面积

东部边界：马山断裂是研究区内规模较大的NW向断层，两盘地层相对错动较大，形成了2个断块地质体，两断块在新周庄以南具隔水性质，分割为相互无水力联系的不同水文地质单元；断层在新周庄至老屯段两盘皆为灰岩，断层为弱透水性质，认为具备作为水文地质边界的条件；断层在老屯以北至前隆段为透水性质，但其影响是局部且有限的。因此，把马山断层作为长清-孝里铺水文地质单元与东部济南泉域水文地质单元的边界。

南部边界：本区的南部和东南部是地表分水岭。在南部，地表分水岭西起黄山岩脉，经贤子峪，冷饭店、兴隆镇南，玉皇顶、木阁寨北等制高点至双山东南出研究区。岩溶水南部补给区以岗辛起向东沿黄石崖、桃花峪分水岭为界。根据水文地质调查证实用地表分水岭近似代替地下分水岭是合理的，构成长清-孝里铺水文地质单元南部的地下水分水岭边界。

西部边界：牛角店断裂为一NW向推测断裂，东南端起自平阴县安城乡大官庄，向东北经栾湾后被第四系覆盖，继续向北延伸。已有地质勘探表明牛角店断裂(黄河以北部分)为透水断裂，故将牛角店断裂(黄河以北部分)作为本区西部的透水边界；黄山岩脉(主支)及其分支阻水性质极为明显，为一阻水边界。基于以上分析故将牛角店断裂(黄河以北段)及黄山岩脉作为长清-孝里铺水文地质单元的西边界。

北部边界：东阿断裂为长清-孝里铺水文地质单元北部边界，西自东阿县高集镇，向东北经齐河县赵官镇、胡官屯镇东部出境。长清-孝里铺水文地质单元边界条件统计情况见表2。

表2 长清-孝里铺水文地质单元边界条件统计

1.2 水文地质背景

1.2.1 水位动态特征

基于2007—2013年曹楼村(CK1)和西关村(CK15)的钻孔地下水水位多年动态变化曲线，可以发现，地下水水位对降水量很敏感，并随降水量的“多—少”呈现“高—低”变化(图2)。根据图2可将区内岩溶地下水水位多年动态变化可分为3个阶段：①2007年1月—2009年12月，水位先下降后上升稳定阶段；②2010年1月—2011年7月，水位先下降后上升再逐渐下降阶段；③2011年8月—2013年12月，水位相对稳定阶段。同时发现地下水水位在枯水期水位总体下降，在丰水期水位总体有所回升。另外CK1、CX15两钻孔点自2007年1月至2013年12月，水位分别上升了0.19m和0.64m，水位呈上升趋势。

1—降水量；2—CK1水位标高；3—CX15水位标高图2 长清-孝里铺水文地质单元岩溶水水位及降水量变化关系曲线图

1.2.2 岩溶发育特征及分布规律

(1)岩溶发育具有层控性

碳酸盐岩是研究区内主要的岩石，具有分布广、厚度大的特点，但各层组碳酸盐岩的生成时代，沉积环境不同，它们的化学成分，矿物成分及结构构造均有一定的差别，故其溶蚀特点及岩溶分布发育程度和特点不尽相同。岩溶最易发育的层位是寒武-奥陶纪三山子组和奥陶纪马家沟群二段、四段和六段和寒武纪中统张夏组，岩溶是易发育的岩性是泥晶灰岩和豹斑灰岩，其次是鲕状灰岩类，第三为白云岩类。由于岩性与层位密切相关，岩溶发育具有一定的层状特点。例如张夏组基本上是巨厚层鲕状灰岩组成，成为单独的岩溶发育的层状溶蚀孔洞-溶隙网络系统，而寒武-奥陶纪三山子组和奥陶纪马家沟群，岩性组成比较复杂多样，虽在岩溶发育上相互联系，但各层均有自己的特点，如其中泥晶和豹斑灰岩以相对均匀性略差的宽溶隙系统为主，白云岩类及角砾灰岩类则以较均匀的层状溶孔及小型孔洞为主，其间又间隔一些岩溶不甚发育的层位。

(2)岩溶发育的分布特征

①水平分布特征。水平分布最主要的特征是由补给区到径流排泄区，岩溶发育呈现由弱到强的变化规律。补给区岩溶发育差，为弱岩溶发育带；补给径流区岩溶发育，为中等岩溶发育带；在径流排泄区，由于地下径流长期作用、构造发育，岩溶十分发育，往往能形成具有开采价值的水源地。

②垂直分布特征。前面已经提到不同岩层的生成时代，沉积环境不同，它们的化学成分，矿物成分及结构构造均有一定的差别，故将岩溶的垂直分布特征分述如下：

奥陶纪马家沟群裂隙、岩溶垂直分布特征：长孝水文地质单元发育深度一般在30～340m之间，其中长孝水文地质单元140～230m这一深度区间发育最普遍，岩溶形态多为溶孔和溶蚀裂隙(图3)。

图3 长清-孝里铺水文地质单元奥陶纪马家沟群裂隙-岩溶发育深度图

寒武-奥陶纪三山子组裂隙-岩溶垂直分布特征：长孝水文地质单元裂隙、岩溶的集中发育深度一般在30～100m和130～200m两个区段，尤其在50～150m这一深度区间，发育最普遍。

寒武纪中统张夏组裂隙-岩溶垂直分布特征：长孝水文地质单元裂隙、岩溶的集中发育深度一般在0～50m和270～320m两个区段。

2 水文地质模型概化

根据地质资料，在考虑最大程度拟合实际水文地质情况下，建立了长清-孝里铺水文地质单元的概念模型及结构模型，利用Visual MODFLOW软件对模型进行数值计算，并根据观测资料对模型进行了为期3年的识别与验证，保证了模型的有效性，随即展开了不同开采条件下的水位预测研究。

2.1 模拟范围

基于前面已详细描述的长清-孝里铺水文地质单元的边界条件，结合前人的研究成果，表明长清-孝里铺水文地质单元基本为一具有系统补给、径流、排泄条件的水文地质单元。故将模拟区总面积确定为780km2，占整个地质单元的83%(图4)。

1—模拟范围；2—透水断层；3—弱透水断层；4—阻水断层；5—岩脉；6—一、二层与三层分界线；7—断层及推测断层界线图4 长清-孝里铺水文地质单元模拟范围示意图

2.2 概念模型

水文地质资料表明研究区主要有4个含水岩组：松散岩类孔隙含水岩组、碳酸盐岩类裂隙岩溶含水岩组、碎屑岩夹碳酸盐岩类岩溶裂隙含水岩组、块状岩类裂隙含水岩组。不同含水岩组的物理性质存在差异，根据研究的重点及目的，将模拟区概化为三层：第一层为孔隙水含水层，为潜水，含水层总厚度10～40m；第二层弱透水层，最大厚度约200m；第三层为岩溶水含水层，具承压性，最大深度为500m左右。

2.3 边界条件概化

2.3.1 孔隙水含水层边界概化

根据前面的地质边界条件可将东部边界南段设为隔水边界，南部、西部、北段边界设为流量边界。

2.3.2 岩溶水含水层边界条件概化

东部边界以马山断裂为界，南端为阻水边界，北段为流量边界。西部边界以西部北端的牛角店断裂为透水边界，黄山岩脉为隔水边界。南部以地下水(地表)分水岭为边界，即零通量边界，地下分水岭与地表分水岭基本一致。北部边界以奥陶纪灰岩顶板埋深600m等值线为界，推测岩溶不发育，可视为阻水边界。

2.3.3 垂向边界

孔隙水含水层自由水面为系统的上边界。通过该边界，潜水与系统外界发生垂向水量交换。孔隙水含水层和岩溶水含水层通过越流进行水量交换，越流量由浅、深层的水头差、弱透水层在垂向上的渗透系数和厚度决定。岩溶水底板作为系统的下边界，底板以下岩层为不透水地层，定义为隔水边界。

3 模型的建立及求解

根据对勘探区水文地质条件的概化，建立数学模型，采用有限差分法进行不同开采条件下的岩溶水开采方案计算，并进行地下水位预报。

3.1 数学模型

地下水流系统具有非均质性、各向同性、空间三维结构、非稳定性，可用如下方程的定解问题来描述：

式中：Kx，Ky和Kz为x，y，z方向渗透系数(m/d)；h为水位标高(m)；W为含水层的源汇项(1/d)；S为含水介质的储水率(1/m)；h0为初始水位(m)；Γ2为二类边界；Kn为边界面法向方向的渗透系数(m/d)；n为Γ2边界的外法线方向；q(x,y,z,t)为二类边界上已知流量函数，流入为正，流出为负，隔水边界为0(m/d)；Ω为渗流区域。

3.2 数值模拟模型

根据模拟区的含水层结构特征、边界条件和地下水流场等，在垂向上分为三层，其中第一层为潜水含水层，二、三层为承压含水层。在平面上剖分成大小为250m×250m的正方形单元，各层均剖分为172行，152列，其中一层有效单元格为6878个，二层有效单元格为6932个，三层有效单元格为12480个。网格剖分示意图见图5。

图5 含水层模型剖分示意图(左：孔隙水；右：岩溶水)

3.3 源汇项的处理

模拟区源汇项有点、线、面3种要素。点状要素主要代表水厂、自备井、工农业开采等源汇项；线状要素主要代表河流(通过River模块处理)、引水渠、灌渠等补给项；面状要素由Recharge模块给出，主要代表降雨入渗、灌溉入渗等补给项。其余源汇项资料通过Well模块带入模型计算，均处理成开采或者补给井。

4 模型的识别与验证

为保证已建立模型的有效性，模型的识别与验证是必不可少的。本研究进行了充分的模型识别与验证。识别期(2010.5.31—2011.6.11)；验证期(2011.6.2—2013.5.26)。模型的识别与检验过程是整个模拟中极为重要的一步工作，本文采用试估—校正法进行了反复修改参数和调整某些源汇项达到了较为理想的拟合结果。经过多次调整参数拟合，得到了3层含水层的水文地质参数。最终的模拟流场与实际流场总体流动方向相同，等值线密间距接近，初步的拟合效果较好。

考虑模拟区水位观测点的分布情况及观测资料的连续性，在第一含水层选择了5个观测点、第三含水层选择了10个观测点，点位具体分布情况见图6。进行水位动态过程曲线拟合，计算各拟合点的计算水位与观测水位的绝对误差和相对误差，并对其进行综合分析，部分观测孔过程线拟合情况见图7，观测值和模拟值的变化趋势基本一致，2条曲线的偏离程度较低，曲线真实的反映了丰枯水期地下水水位的动态变化，基本反映了地下水系统的水力特征，模型的有效性得到了验证。

1—岩脉；2—透水断层；3—弱透水断层；4—阻水断层；5—断层及推测断层界线；6—孔隙水观测井；7—岩溶水观测井图6 地下水位监测孔分布示意图

5 开采方案模拟与预测分析

5.1 开采方案模拟

在已建立的模型上，进行不同开采方案的开采动态响应研究，开采的目的层位为岩溶含水层，在研究区内选择了一处富水地段，根据原有地质资料显示的岩溶水储量情况确定了2种开采方案：

表3 开采方案具体实施方案

在进行多方案预测之前，进行了自然均衡状态下曹楼抽水试验和相同地点的少量增采试验，对两者的地下水流场分析，证明了长孝水文地质单元有一定增采潜力。在将每年分为枯水期与丰水期2个应力期的基础上，分别预测分析了2个方案运行10年后地下水流场变化情况。预测期间，2个方案降水量均采用1995—2004年枯水系列年资料，模型中以559m平均降雨量进行预测。

5.2 开采方案对地下水位的影响

运行方案1后，地下水位产生明显变化，在布井地段出现了一个明显的地下水降落漏斗，在曹楼附近出现最大降深约9.02m，年降幅为0.902m；运行方案2后，地下水位降落漏斗进一步扩大，曹楼最大降深为12.25m，年降幅达1.225(图8)。随着开采量增加，岩溶水得不到有效补给，导致降落漏斗增大，水位降深加大。运行两方案后，由模拟区马山断裂长清北排泄区观测孔水位变化情况可以看出，观测孔水位有持续下降趋势，虽然下降幅度不大，但改变了北边界的补排方向。增加的开采方案不仅增加了地下水补给量，同时也减少了排泄量减少，但补给增量是有限度的，故导致了地下水水位的降低。

1—观测值；2—模拟值图7 观测孔水位过程线拟合图(左：识别期；右：验证期)

图8 岩溶水预测流场示意图(左：方案1；右：方案2)

通过数值法对岩溶水开采方案模拟分析表明，实施开采方案后，地下水水位变化明显，在开采井附近产生地下水降落漏斗，对水文地质单元的地下水位也有一定影响。枯水系列年的2种开采方案改变了地下水均衡状态。自然情况下，马山断裂为岩溶水的排泄边界，而实施开采方案后，激化增加补给量，袭夺部分边界排泄量，但激化与袭夺总量是有限度的，结合该区的实际情况，各方案说明在连续枯水系列年条件下，长清地区开采量不宜超过14万m3/d(含分散开采)。若考虑东边界以东外围水源地的影响，应建立长孝-济南数值模型，统一进行地下水资源分析与预测，才能进行精度更高的模型研究。