济南市历城区地下水数值模拟与预测

朱恒华，张晓伟，王 玮，贾 超，刘治政，刘柏含

(1.山东省地质调查院，山东 济南 250013；2.中国地质大学(武汉)，湖北 武汉 430074；3.山东大学 土建与水利学院，山东 济南 250061； 4.山东大学 海洋研究院，山东 青岛 266237)

随着社会经济的不断发展，城市工业和农业用水量日益增加，水资源量消耗剧增。由于地表水易受污染，利用率较低，开采地下水的需求日益凸显。然而，过度开采地下水会造成地面沉降、岩溶塌陷等地质问题[1]。因此，科学合理地开采地下水就显得尤为重要。近年来GMS被广泛地应用于地下水资源的计算评估预测及相关的环境水文地质研究中并取得了良好的效果[2，3]。在地下水数值模型建模过程中，相关水文地质参数对模型精度影响较大，其中又以渗透系数最为敏感[4]。PEST是GMS中自动调参模块，因此选取渗透系数进行参数自动估计。运用PEST来调整选定的参数，重复运行MODFLOW模块的计算，直到计算结果和野外观测值相吻合。历城区位于济南市，前人研究多集中于岩溶水，对松散岩类孔隙水研究较少。本文通过GMS软件建立历城区浅层地下水数值模型，对不同开采方案下的地下水流场进行预报，为浅层地下水资源开采提供参考依据。

1 研究区地质概况

历城区位于济南市中部，地势南高北低，如图1。根据区内地貌特征，可划分为堆积地貌和侵蚀地貌，本区南部是一个以古生代地层为主体的北倾单斜构造，北部则是沉积了厚度较大的第四系地层和第三系松散堆积物的凹陷，区内地层由老到新依次为古生界寒武系、奥陶系、石炭-二叠系；新生界第三系及第四系。

图1 研究区地理位置图

研究区处于鲁中山地的北缘以及黄河下游冲洪积平原区，区内地层、构造条件复杂，根据地层岩性、含水介质及地下水运动、储存等特征，区内地下水类型可分为松散岩类孔隙水、碎屑岩类孔隙裂隙水、碳酸盐岩类裂隙岩溶水。除基岩裸露区外，松散岩类孔隙水在全区普遍分布。含水层岩性为粉土、粉细砂、中细砂、粗砂、砾石及粘质砂土等，地下水多属潜水或微承压水，如图2；碎屑岩类孔隙裂隙水主要分布于北部、西北部，孔隙裂隙不甚发育，富水性较差；碳酸盐岩类裂隙岩溶水主要赋存在奥陶系灰岩中，分布于南部及黄河以北局部地段，与浅层孔隙水水力联系微弱。

图2 黄河冲积平原区水文地质剖面图

2 地下水流模型

2.1 水文地质概念模型

2.1.1 地层概化

研究区主要含水层确定为浅层孔隙水含水层，第四系厚度20～170 m，总体呈现由南向北渐厚的趋势。浅层孔隙水以下为中层咸水，不进行地下水开采，故将含水层概化为深度170 m的非均质、各向同性的潜水微承压含水层。

2.1.2 边界条件

1)垂向边界：上边界接受大气降水、排泄方式为大气蒸发，是不断变化的水量交换边界，定为流量边界；下边界因水量交换微弱，概化为相对隔水边界。

2)侧向边界：研究区南部为基岩出露区，概化为隔水边界[5]；北部以李家断层为界概化为隔水边界，其余边界作为流量边界，根据达西定律计算截面流量。

2.1.3 源汇项

(1)大气降水补给量以面状补给赋值，根据降雨入渗系数法[6]计算本区松散岩类孔隙水大气降水补给总量为4 249.965万 m3/a(合11.64万 m3/d)。灌溉回渗量以面状补给赋值，依据研究区内的灌溉定额、类型及包气带岩性等计算区内灌溉回渗补给总量为187.663万 m3/a，其中井灌回渗补给量为59.251万 m3/a，黄灌回渗补给量为132.973万 m3/a；区内对地下水具有一定补排规模的河渠有黄河、邢家度总干渠。河渠渗漏量以线状补给赋值[7]。根据测流数据，黄河断面的单侧单宽补排量为1.102 m3/d·m，侧渗补给总量为172.009万 m3/a；邢家渡干渠渗漏系数为1.010 m2/d·m，侧渗补给总量为52.803万 m3/a。

(2)根据区内水位埋深和包气带岩性计算本区潜水蒸发总量为3 650.671万 m3/a；区内地下水开采量为775.818万 m3/a；区内主要排泄河流为小清河，经计算，本区地下水排泄小清河总量为48.681万 m3/a。地下水侧向径流排泄总量为137.547万 m3/a。

2.2 数学模型

对于非均质各向同性的三维非稳定流地下水模型，其数学方程描述为：

(1)

定解条件：

H(x，y，z)|t=0=H0(x，y，z)，(x，y，z)∈Ω

(2)

H(x，y，z，t)|(x，y，z)∈B1=H1(x，y，z，t)，(x，y，z)∈B1，t>0

(3)

(4)

式中：Kxx，Kyy，Kzz为xyz方向的渗透系数，LT-1；h为水头，L；w为单位体积含水层在单位时间流出或流入地下水的体积，T-1；Ss为储水系数；Ω为初始条件，即给定边界所有点的定水头；B1为第一类边界条件；B2为第二类边界条件。

3 地下水数值模型

3.1 模型的建立

采用GMS中MODFLOW模块建立历城区地下水流数值模型，将研究区空间离散为150×150个单元格，其中有效单元格19 726个。根据已有资料确定2018年6月-2020年6月为模拟期，时间步长为30 d，共24个应力期。导入高程数据，输入各源汇项数据。研究区网格剖分如图3。

图3 研究区三维网格剖分图

根据研究区的地质构造、岩性及富水性将研究区按水文地质参数进行分区[8,9]，如图4，初始渗透系数和给水度取值见表1。

图4 研究区参数分区图

表1 水文地质参数初值表

根据研究区内观测水位进行插值拟合出2018年6月浅层孔隙水流场作为初始流场，如图5。

图5 研究区初始流场拟合图

3.2 模型识别检验

传统的水文地质参数分区存在较强的主观性，忽略了渗透系数的空间差异性以达到局部最优解，而PEST模块在初始参数分区的基础上可通过添加散点反演渗透系数场从而达到全局最优解。

选取2019年1月至2019年12月观测井的水位数据进行PEST参数反演，所得渗透系数场如图6，观测井空间分布如图7，并以此进行地下水动态检验。典型长观井水位的实测值与计算值拟合情况如图8，可以看出实测值与模拟计算值动态变化趋势一致，拟合情况较好，可以此模型进行地下水动态预测[10]。

图6 渗透系数反演场 图7 研究区观测井空间分布图

图8 研究区典型长观井水位过程线拟合图

4 地下水动态及流场预测

地下水开采应以可持续发展的观点，在兼顾社会经济和环境保护的基础上，合理有效的开发利用水资源。因此，在已有资料的基础上，设计不同的开采方案，通过运行地下水数值模型来预测不同开采方案下的地下水流场。

4.1 方案设计

从充分利用地下水资源的角度考虑，在现状开采量2.13万 m3/d的基础上分别增采10%，20%，50%，70%，预测5 a后的浅层地下水流场。预测起始时间为2020年6月，目标时间为2025年6月。预测期间，大气降水按多年平均降水量647.9 mm取值。

4.2 结果分析

(1)经过模型预测，在现状开采量的基础上增采10%的情况下，5年后浅层地下水流场见图9，研究区内典型长观孔水位动态变化见图13(a)。可以看出，当设计开采量为2.343万 m3/d时，地下水流场水位变幅在0.1～0.4 m之内，地下水降落漏斗区域面积减小，水位有明显抬升。

图9 增采10%地下水流场

(2)在现状开采量的基础上增采20%的情况下，5年后浅层地下水流场见图10，研究区内典型长观孔水位动态变化见图13(b)。相比于增采10%，研究区内典型长观孔水位抬升变缓，水位变幅在0.1～0.4 m之间。地下水降落漏斗区域面积仍在减小。

图10 增采20%地下水流场

(3)在现状开采量的基础上增采50%的情况下，5年后浅层地下水流场见图11，研究区内典型长观孔水位动态变化见图13(c)。当设计开采量达到3.195万 m3/d时，研究区内典型长观孔水位于开采4 a后出现水位下降的趋势，地下水流场水位变幅处于0.05～0.35 m之间。

图11 增采50%地下水流场

(4)在现状开采量的基础上增采70%的情况下，5年后浅层地下水流场见图12，研究区内典型长观孔水位动态变化见图13(d)。可以看出，当设计增采量达到1.49万 m3/d时，研究区内典型长观孔水位于开采2 a后出现水位下降的情况，且下降趋势较增采50%更为明显。地下水流场水位变幅仍处于0.05～0.35 m之间。地下水降落漏斗中心区域开始出现水位下降的情况，可以预见，若持续以设计70%增采量开采浅层地下水将会出现降落漏斗面积扩大，研究区水位持续下降等环境问题。

图12 增采70%地下水流场

图13 不同开采方案下典型长观孔未来5a地下水动态变化曲线

5 结语

通过数值模型预测，将历城区浅层孔隙水开采量控制在3.195万 m3/d之内，可以有效地控制研究区浅层地下水水位下降，防止地下水降落漏斗持续扩大，进而引起地面沉降和地下水污染等环境地质问题，达到地下水资源充分且可持续利用的目标。