基于Visual MODFLOW地下水数值模拟的水源地保护研究
曹振东1,3, 谭廷静1,3, 杨明星2, 宋小庆1,3, 蒲秀超1,3
1.贵州省地质矿产勘查开发局111地质大队,贵州 贵阳 550081
2.贵州理工学院资源与环境工程学院,贵州 贵阳 550003
3.贵州地质工程勘察设计研究院有限公司,贵州 贵阳 550081
通信作者简介: 杨明星(1986—),男,副教授,主要从事岩溶水文地质方面的研究工作。Email: yangmingxing@git.edu.cn

第一作者简介: 曹振东(1984—),男,高级工程师,主要从事地下水资源保护方面的研究工作。Email: caozd16@mails.jlu.edu.cn

摘要

明确地下水运移条件是城市地下水水源地保护的基础。为确保沈阳市某地下水源地的可持续开发,依据沈阳市某地下水水源地的水文地质条件,构建了研究区水文地质概念模型和地下水数值模型,采用Visual MODFLOW软件对模型进行求解,并基于大量监测井水位数据对模型进行参数的识别和验证。依据校正后的数值模型分析了水源地运营后地下水漏斗的范围及水位降深,并对水源地停采后水位恢复状况进行了模拟预报。结果表明: 通过模型计算可知,研究区水源地地下水含水层总补给量为62 230 m3/d,总排泄量为63 400 m3/d,均衡差为-1 170 m3/d,多年地下水动态呈负均衡状态; 通过对研究区水源地开采2 a、5 a、10 a后承压水状态进行预测,地下水水位分别平均下降6 m、8 m、9 m,形成中心漏斗区,面积分别为54.56 km2、65.04 km2、65.80 km2,开采初期降落漏斗急剧扩张,然后速度逐渐放缓; 研究区水源地承压含水层开采对周围流场产生了一定的影响,但这种影响主要在开采期较为明显,停抽1 a后漏斗区逐渐恢复,为了水源地的可持续利用及保护,建议严格控制开采量,加强地下水的监测与管理。预测结果可以为水源地保护提供技术支撑。

关键词: 地下水; 水源地保护; 数值模拟; Visual MODFLOW
中图分类号:P641.2 文献标志码:A 文章编号:2095-8706(2023)05-0091-08
Research on water source protection based on Visual MODFLOW groundwater numerical simulation
CAO Zhendong1,3, TAN Tingjing1,3, YANG Mingxing2, SONG Xiaoqing1,3, PU Xiuchao1,3
1. No.111 Geological Party, Guizhou Bureau of Geology and Mineral Exploration & Development, Guizhou Guiyang 550081, China
2. School of Resource and Environment Engineering, Guizhou Institute of Technology, Guizhou Guiyang 550003, China
3. Guizhou Geological Engineering Survey, Design and Research Institute Co., Ltd., Guizhou Guiyang 550081, China
Abstract

Clarifying the migration conditions of groundwater is the key basis for the protection of urban groundwater sources. In order to serve the sustainable development of the groundwater source in Shenyang, the authors constructed the hydrogeological conceptual model and groundwater numerical model according to the hydrogeological conditions of the study area. The Visual MODFLOW software was used to solve the model, and the model parameters were identified and verified based on the water level data of observation well. The range of groundwater funnel and the depth of water level drop after water source operation were analyzed according to the corrected numerical model, and the water level recovery after the stop-mining of water source was simulated and predicted. The results showed that the total recharge of groundwater aquifer in the study area was 62 230 m3/d, and the total discharge was 63 400 m3/d, with equilibrium difference of -1 170 m3/d. And groundwater dynamic was in a negative equilibrium state for many years. Through the prediction of the confined aquifer situation after 2 a, 5 a and 10 a of water source exploitation, the water level decreased by an average of 6 m, 8 m and 9 m respectively, and the area of the central funnel area was 54.56 km2, 65.04 km2 and 65.80 km2 respectively. The funnel expanded rapidly in the early stage of mining, and then the speed gradually slowed down. The exploitation of confined aquifer in the water source had a certain impact on the surrounding flow field, but this impact was mainly obvious during the mining period, and the funnel area gradually recovered after one year of stop-pumping. It is suggested that the exploitation amount should be strictly controlled and the monitoring and management of groundwater should be strengthened for the sustainable utilization and protection of water source. The prediction results could provide some technical support for water source protection.

Keyword: groundwater; water source protection; numerical simulation; Visual MODFLOW
0 引言

作为城市供水的重要水源, 地下水是区域民生、经济安全发展的基础保障。从比例上分析, 全国水资源总量为29.6×1011 m3。其中, 地下水资源量为8.20×1011 m3, 占我国水资源总量的27.7%[1]。因含水介质类型、埋藏条件等赋存条件差异, 导致地下水补给、径流、排泄过程较为复杂, 难以准确刻画, 是制约地下水水源地保护的关键因素[2, 3]。针对此项工作, 国内外采用了多种手段对地下水赋存条件进行分析。其中, 在建立准确地下水概念模型及设定初始条件的基础上, 利用计算机数值模拟计算技术, 可实现对地下水运移过程的模拟[4, 5]。Visual MODFLOW软件由加拿大Waterloo水文地质公司开发, 由于特定的模块化程序结构、简单化的离散方法以及多样化的求解方法等优点, 得到广泛的应用[6]。韩琳等[7]和吴欣等[8]应用Visual MODFLOW对地下水进行了数值模拟研究及水位预测; 马佳[9]利用Visual MODFLOW模拟了唐山市乐亭县工业园区地下水的溶质运移情况, 为工业园区地下水开发利用做出了参考; 刘天宇[10]则对新乡市的某水源地进行了数值模拟与分析, 进行了合理的水资源评价。此外, 地下水数值模拟在矿山排水、地铁抗浮水位和地面沉降等方面也得到了广泛应用[11, 12, 13, 14]。本文在沈阳市某水源地基础地质条件分析的基础上, 根据研究区2018 年41 口水井的动态资料(水位、流量), 利用Visual MODFLOW地下水数值模拟软件构建该水源地数值模型, 通过数学模型刻画、参数率定、参数的识别验证等过程, 建立了研究区地下水数值模型, 并模拟了不同运营情景下地下水水流过程。在此基础上, 得到了开采后降落漏斗的发展情况, 提出了地下水开采管理建议, 旨在为水源地保护工作提供了有力的技术支撑。

1 研究区概况

研究区位于沈阳经济技术开发区内, 所在地区地势平坦, 地面高程为33.5~19.5 m, 由北东向南西微倾, 位于浑河新老冲洪积扇的中部, 地貌由北向南可分为一级阶地、高河漫滩、低漫滩。开发区境内流经的河流主要是细河, 全长78.2 km, 位于沈阳市西南部城市出口处, 源于铁西区卫工河南端, 流经铁西区、于洪区、开发区、辽中区, 在辽中区黄腊坨子村汇入浑河(图1)。

图1 研究区位置Fig.1 Location of the study area

研究区属于半湿润暖温带气候型, 受西北欧亚大陆和东南海洋性高压气团控制, 四季分明, 气候宜人。年平均气温7~8 ℃, 多年平均降水量为708.6 mm。

研究区位于阴山EW向复杂构造带东延部位, 主要被大厚度第四系松散堆积物掩埋。EW向构造、新华夏系构造均为压性断裂, 华夏系构造有压性断裂和褶皱。

2 地下水数值模型的建立
2.1 水文地质概念模型建立

研究区水源地有水厂3处, 水源井41眼: 一水厂水源井33眼; 二水厂水源井4眼, 位于二水厂院内及灌渠路; 三水厂拟建水源井4眼, 位于沈阳市经济技术开发区三水厂内、开发北六路沿线和四环路沿线一带。水源井运行后日取水量为1万m3/d, 为小型承压水水源地。本次研究范围南部以浑河为界, 北侧到大王庄村—沙岭村一带, 西侧到马贝村—二牤牛村, 东侧至北李官村—郑家村一带, 研究区总面积约154.56 km2

根据区域水文地质资料的调查分析[15], 研究区含水层自上而下分别为全新统松散岩类孔隙潜水、上更新统松散岩类孔隙承压水和下更新统承压水, 两者之间为稳定的隔水层, 岩性为粉质黏土, 厚度较大。地面以下20 m左右有厚度约10 m的粉质黏土层, 起到隔水作用(粉质黏土层厚1.70~17.6 m, 最厚达26.5 m, 顶板埋深为20~30 m), 实现了全新统松散岩类孔隙潜水与上更新统松散岩类孔隙承压水的隔离。此外, 第四系中更新统以冲洪积为主, 次为冲积成因, 岩性以黄褐色、灰黄色粉质黏土为主, 厚度为10~22 m, 是上更新统松松散岩类孔隙承压水的隔水底板。其中, 上更新统松散岩类孔隙承压水为水源井的主要开采层位, 含水层岩性为中粗砂、砾砂, 少量为圆砾层, 含水层厚25.02~48.60 m。承压水位埋深为4~10m, 静水位年变幅为0.8~2.5 m。抽水试验结果表明, 上更新统含水层渗透系数为30~50 m/d, 单井出水量为3 000~3 500 m3/d。

2.1.1 边界条件及水力特征概化

垂向边界的概化。依据基础水文地质条件分析, 设定研究区承压含水层隔水顶板为研究区的上边界, 通过该边界, 上部潜水含水层与承压含水层发生垂向越流; 研究区的下部为承压水的隔水底板, 厚度较大, 岩性为粉质黏土, 概化为隔水边界。

侧向边界的概化。依据研究区地下水天然流场特征, 地下水由北东向南西径流, 北东及南西侧概化为流量边界, 分别为承压水的流入与流出边界; 天然条件下研究区西北侧边界方向与地下水流线方向接近, 概化为零流量边界; 东南侧为浑河, 虽然河床未切割承压含水层, 但承压水位在河道处水位抬升为水丘, 可概化为隔水边界。水源井运行后, 预报模型中边界条件将依据开采流场进行调整。

研究区地下水系统符合质量守恒定律, 含水层分布广, 在常温常压下地下水运动符合达西定律。研究区范围相对较小, 水文地质参数随空间变化较小, 且没有明显的方向性, 所以可以将研究区水文地质参数概化成均质各向同性。由于地下水系统渗流运动要素随时空变化, 故地下水含水系统概化为非稳定流。

综上所述, 研究区地下水系统可概化为均质、各向同性、二维非稳定地下水流系统。

2.1.2 参数分区

研究区为均质含水层, 渗透系数通过抽水试验的结果获取, 区内取值30 m/d。通过查阅该水源地供水水文地质调查报告, 贮水系数取值2.6×10-3

2.2 地下水数值模型建立

2.2.1 数学模型建立

根据研究区的均质、各向同性及空间二维结构的特征及非稳定承压水系统的划分, 数值模拟将采用如下微分方程(1)的定解问题来描述

xThx+yThy+ε=μ·ht (x, y)D, t0h(x, y, t)|t=0=h0(x, y) (x, y, z)DKhhn|Γ1=q(x, y, t) (x, y, z)Γ1, t0 。(1)

式中: D为地下水系统的模拟渗流区; h为地下水水位, m; T为含水层导水系数, m2/d; K为渗透系数, m/d; μ为贮水系数; t为时间, d; Γ 1为给定流量边界; n为边界的外法线方向; h0(x, y)为初始地下水位, m; q(x, y, t)为给定流量边界上的单宽流量, m2/d; ε为垂向越流量, m3/d。

2.2.2 网格剖分

本次模拟的模拟面积为154.56 km2。根据Visual MODFLOW的要求, 在水平方向上用相互平行的铅垂格线将研究区承压含水层剖分成86行、91列的柱形网格。单元格顶面为上隔水层与含水层的边界, 单元格底面为下隔水层与含水层的边界, 网格单元横截面尺度约200 m×200 m, 单元格总数为7 826个, 其中有效单元格4 022个。

2.2.3 模型识别与验证

模型识别和验证是数值模拟极为重要的过程, 通常需要进行多次的参数调整与运算。运行模拟程序, 可得到概化后的水文地质概念模型, 从而模拟出给定水文地质参数和各均衡项条件下的地下水流场空间分布。通过拟合同时期的实际地下水流场, 识别和验证水文地质参数、边界值和其它均衡项, 使建立的模型更加符合研究区的实际水文地质条件。以2018年7月1日到2018年12月31日作为模型的识别验证期(共计184 d), 由地下水实际观测资料提取研究区实际地下水流场数据, 计算建模所需源汇项, 按照软件操作流程输入各含水层参数, 最终通过模拟计算结果与实际观测结果的对比, 确定符合实际的研究区概化水文地质数值模型。

以2018年7月1日承压水实际观测水位为模型识别验证过程的初始水位, 研究区承压含水层初始流场见图2。

图2 研究区识别验证期承压含水层初始流场Fig.2 Initial flow field diagram for confined water during the identification and verification period of the study area

由于参数初值的选取客观地反映了实际的水文地质条件, 加之细致地调参拟合, 模型识别验证取得了较好的结果。本文共选取了6口同一承压水层的观测井作为识别和检验的校准数据点及水位拟合情况(图3)。以观测井G1为例展示拟合的情况(图4), 由于观测井本身为开采井, 观测水位反映了地下水开采的影响, 因此, 该模拟精度能够满足生产需求, 该模型反应了实际情况下水位随时间的变化特征。由图4可见, 观测井G1在各个时段的模拟水位与观测水位拟合程度较好, 水位拟合差小于0.65 m, 识别结果可靠, 所建立的模拟模型基本反映了研究区地下水系统的水力特征, 达到了模型精度要求。识别验证期各源汇项均衡结果见表1

图3 研究区识别验证期承压水观测井位置分布(左)及承压水水位拟合图(右)Fig.3 Locations distribution of observation wells of confined water (left) and confined water level fitting diagram (right)during the identification and verification period

图4 G1观测井承压水水位拟合Fig.4 Water level fitting diagram for G1 observation well of confined water in the study area

表1 研究区识别验证期各均衡项计算结果(单位: m3/d) Tab.1 Calculation results of each equilibrium term during the identification and verification period of the study area (unit: m3/d)
2.3 地下水水流预测

以识别和校正后的模型为基础, 源汇项的处理方法与识别期相同, 其中侧向补排量以水源地运行的流场状况进行计算, 水源井的开采量按设计要求加载到模型中[11, 12]。预测模型的初始水位取2019年1月1日统测水位, 预测时段根据水源地服务年限进行预报。

依据水源地规划, 水源井的服务年限为10 a, 即2019—2028年。由于承压含水层上部的潜水含水层直接接受大气降水补给, 地下水补径排速率较快, 且潜水地下水开采较少, 因此水源井开采过程中, 虽开采层接受上部潜水含水层越流补给, 但对地下水位影响有限。因此, 在预报期仅选择水源地开采2 a、5 a、10 a及停采后, 开采层即承压含水层地下水流场及水位变化进行分析。

通过2 a后(2020 年)的承压含水层流场(图5)可以看出, 水源井运行后, 研究区承压水水位在9.6~22.2 m之间, 在水源井周边形成明显的区域降落漏斗, 地下水流向也发生了明显变化, 由原来的自北东向南西流向, 变化为四周向漏斗中心汇流。与初始水位相比, 水源地运行后研究区内承压水水位平均下降6 m, 其中: 一水厂水源井周边水位下降8~11 m, 二水厂水源井周边承压水水位下降6~7.1 m, 三水厂水源井周边承压水水位下降3~4.3 m。区内承压水最大水位下降幅度为11.2 m, 位于一水厂。

图5 研究区水源地运营2 a后承压含水层流场(左)及降深等值线(右)Fig.5 Confined aquifer flow field (left) and drawdown contour (right) after 2 years of water source operation

通过5 a 后(2023 年)的承压含水层流场(图6)可以看出, 水源井运行后, 研究区水源地内承压水水位在9.2~22.0 m之间, 水源井周边的区域降落漏斗进一步扩大, 整个水源地承压水水位出现一定幅度下降。与初始水位相比, 水源地运营后区内承压水水位平均下降8 m, 最大水位下降12.1 m, 位于一水厂周边。

图6 研究区水源地运营5 a后承压含水层流场(左)及降深等值线(右)Fig.6 Confined aquifer flow field (left) and drawdown contour (right) after 5 years of water source operation

通过10 a后(2028 年)的承压含水层流场(图7)可以看出, 水源井运行后, 研究区水源地内承压水水位在9.2~22.0 m之间, 水源井周边的区域降落漏斗进一步扩大, 整个水源地承压水水位出现一定幅度下降, 但漏斗中心水位趋于稳定。与初始水位相比, 水源地运营后区内承压水水位平均下降9 m, 最大水位下降12.1 m, 位于一水厂周边。10 a后, 依据最大降落漏斗中心(一水厂水源井)的下降速率计算, 年平均下降速率为1.21 m/a。二水厂及三水厂水源井的漏斗中心下降速率相对较慢, 二水厂年平均下降速率为0.72 m/a。三水厂年平均下降速率为0.80 m/a。

图7 研究区水源地运营10 a后承压含水层流场(左)及降深等值线(右)Fig.7 Confined aquifer flow field (left) and drawdown contour (right) after 10 years of water source operation operation

3 水源地保护研究

以降深5 m等值线圈定的范围作为降落漏斗区, 对水源地地下水开采引起的地下水降落漏斗变化情况进行分析。开采2 a后(2020 年), 地下水降落漏斗的面积为54.56 km2, 占研究区总面积的35.3%; 开采5 a后(2023年), 地下水降落漏斗的面积为65.04 km2, 占研究区总面积的42.1%; 开采10 a后(2024年), 地下水降落漏斗的面积为65.80 km2, 占研究区总面积的42.6%。由分析结果可以发现, 开采初期降落漏斗急剧扩张, 然后速度逐渐放缓, 开采10 a后漏斗总面积未达到总模拟区的50%。

由水源地停采1 a后的承压含水层流场(图8)可以发现, 停采后, 研究区承压水水位出现恢复, 降落漏斗迅速消失, 之后水位呈现整体恢复。停采1 a 后研究区承压水水位在16~23 m。

图8 研究区水源地停采1 a后承压含水层流场(左)及降深等值线(右)Fig.9 Confined aquifer flow field (left) and drawdown contour (right) after 1 year of stop-mining in the water source

根据以上地下水数值模拟结果, 研究区水源地开采会对周边地下承压含水层流场产生一定程度的影响, 且随着开采时间延长, 地下水降落漏斗范围逐渐扩大, 漏斗中心水位逐渐趋于稳定。研究区水源地开采后, 水源地范围内相同开采层内无其他农用井, 不会引起周边水井的疏干, 停采后, 水位恢复速度较快, 因此整个开采过程对周边地下水位的影响较小。但是, 仍要按照水源地保护区等级一级去严格管理, 不可持续性超采, 要严格保护、监测、监督管理, 防治水污染, 防范水污染事故发生, 保障饮用水源水量、水质安全。

4 结论

(1)基于研究区地质及水文地质条件, 构建了地下水数值模型, 计算出含水层总补给量为62 230 m3/d, 总排泄量为63 400 m3/d, 均衡差为-1 170 m3/d, 表明研究区水源地地下水动态呈负均衡状态。

(2)通过对研究区水源地开采2 a、5 a、10 a后的承压含水层流场进行预测, 地下水水位分别平均下降6 m、8 m、9 m, 形成中心漏斗区, 面积分别为54.56 km2、65.04 km2、65.80 km2, 开采初期降落漏斗急剧扩张, 然后速度逐渐放缓。

(3)研究区水源地承压含水层开采对周围流场产生了一定的影响, 但这种影响主要在开采期较为明显, 研究区水源地停抽1 a后漏斗区逐渐恢复。为了水源地的可持续利用及保护, 建议严格控制开采量, 并加强对水源地地下水资源的监测与管理。

(责任编辑: 刘丹, 王晗)

参考文献
[1] 中华人民共和国水利部. 2021年中国水资源公报[R]. 北京: 中华人民共和国水利部, 2021.
Ministry of Water Resources of the People's Republic of China. 2021 China Water Resources Bulletin[R]. Beijing: Ministry of Water Resources of the People's Republic of China, 2021. [本文引用:1]
[2] 许颢砾, 王大庆, 邓正栋, . 基岩岛屿地下水流场数值模拟研究[J]. 中国地质调查, 2020, 7(4): 95-103.
Xu H L, Wang D Q, Deng Z D, et al. Research on the numerical simulation of the groundwater flow field in bedrock island s[J]. Geological Survey of China, 2020, 7(4): 95-103. [本文引用:1]
[3] 卢洪健, 王卓然. 地下水模拟方法与应用软件研究进展[J]. 地下水, 2022, 44(6): 49-52.
Lu H J, Wang Z R. Research progress in groundwater simulation methods and application softwares[J]. Ground Water, 2022, 44(6): 49-52. [本文引用:1]
[4] 王浩, 陆垂裕, 秦大庸, . 地下水数值计算与应用研究进展综述[J]. 地学前缘, 2010, 17(6): 1-12.
Wang H, Lu C Y, Qin D Y, et al. Advances in method and application of groundwater numerical simulation[J]. Earth Science Frontiers, 2010, 17(6): 1-12. [本文引用:1]
[5] 王大庆, 许颢砾, 邓正栋, . 基岩岛屿地下水数值模拟发展研究现状[J]. 中国地质调查, 2019, 6(3): 68-74.
Wang D Q, Xu H L, Deng Z D, et al. Development and research status of numerical simulation on the groundwater of bedrock island s[J]. Geological Survey of China, 2019, 6(3): 68-74. [本文引用:1]
[6] Haque M A M, Jahan C S, Mazumder Q H, et al. Hydrogeological condition and assessment of groundwater resource using visual modflow modeling, Rajshahi city aquifer, Bangladesh[J]. Journal of the Geological Society of India, 2012, 79(1): 77-84. [本文引用:1]
[7] 韩琳, 颜翠翠, 赵振伟, . 基于Visual MODFLOW的地下水数值模拟研究及水位预测[J]. 山东国土资源, 2021, 37(3): 67-74.
Han L, Yan C C, Zhao Z W, et al. Groundwater numerical simulation and water level prediction based on Visual MODFLOW[J]. Shand ong Land and Resources, 2021, 37(3): 67-74. [本文引用:1]
[8] 吴鑫, 孙伯明, 陈菁, . 基于Visual MODFLOW的挠力河流域地下水数值模拟与预测分析[J]. 水电能源科学, 2020, 38(12): 37-40, 23.
Wu X, Sun B M, Chen J, et al. Numerical simulation and prediction analysis of groundwater in Naoli river basin based on Visual MODFLOW[J]. Water Resources and Power, 2020, 38(12): 37-40, 23. [本文引用:1]
[9] 马佳. 基于MODFLOW的乐亭县工业聚集区地下水数值模拟研究[D]. 南昌: 东华理工大学, 2020.
Ma J. Numerical Simulation of Groundwater in the Industrial Gathering Area of Laoting County Based on MODFLOW[D]. Nanchang: East China University of Technology, 2020. [本文引用:1]
[10] 刘天宇. 新乡市原武—包厂水源地地下水数值模拟与水资源评价[D]. 北京: 中国地质大学(北京), 2017.
Liu T Y. Numerical Simulation and Water Resources Evaluation of Yuanwu-Baochang Groundwater Source in Xinxiang[D]. Beijing: China University of Geosciences (Beijing), 2017. [本文引用:1]
[11] 李贵仁, 赵珍, 刘大金, . 基于Visual MODFLOW的复杂大水矿山地下水数值模拟[J]. 勘察科学技术, 2019(3): 34-38.
Li G R, Zhao Z, Liu D J, et al. Groundwater numerical simulation of complex and water abundant mine based on Visual MODFLOW[J]. Site Investigation Science and Technology, 2019(3): 34-38. [本文引用:2]
[12] Tiwary R K, Dhakate R, Rao V A, et al. Assessment and prediction of contaminant migration in ground water from chromite waste dump[J]. Environmental Geology, 2005, 48(4-5): 420-429. [本文引用:2]
[13] 南天, 曹文庚, 王卓然, . 利用趋势化随机参数场的地下水流数值模拟优化方法[J]. 现代地质, 2022, 36(2): 591-601.
Nan T, Cao W G, Wang Z R, et al. Optimized groundwater numerical simulation model with trending parameter field[J]. Geoscience, 2022, 36(2): 591-601. [本文引用:1]
[14] Ye S J, Luo Y, Wu J C, et al. Three-dimensional numerical modeling of land subsidence in Shanghai, China[J]. Hydrogeology Journal, 2016, 24(3): 695-709. [本文引用:1]
[15] 王瑾. 沈阳地铁一号线一标段水文地质条件分析[J]. 科技信息, 2013(3): 93-94.
Wang J. Analysis of hydrogeological conditions of the first section of Shenyang Metro Line 1[J]. Science and Technology Information, 2013(3): 93-94. [本文引用:1]