基于Visual MODFLOW地下水数值模拟的水源地保护研究

曹振东, 谭廷静, 杨明星, 宋小庆, 蒲秀超

(1. 贵州省地质矿产勘查开发局111地质大队,贵州 贵阳 550081; 2. 贵州理工学院资源与环境工程学院,贵州 贵阳 550003; 3. 贵州地质工程勘察设计研究院有限公司,贵州 贵阳 550081)

0 引言

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

(1)

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。

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

图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

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后漏斗区逐渐恢复。为了水源地的可持续利用及保护,建议严格控制开采量,并加强对水源地地下水资源的监测与管理。