基于MODFLOW的挠力河流域地下水承载力综合评价

2024-03-08郗鸿峰刘春伟杜新强刘小天马河宽于令芹

地下水 2024年1期

郗鸿峰,刘春伟,杜新强,刘小天,马河宽,于令芹

(1.山东省地质矿产勘查开发局八〇一水文地质工程地质大队,山东济南 250000;2.吉林大学新能源与环境学院,吉林长春 130000)

目前,地下水承载力综合评价在全国各个地区得到了广泛应用,各学者也构建了适合不同区域的评价指标体系及评价等级;在东北农灌区应用较广的8个指标为地下水资源耕地灌溉率、地下水开发利用程度、地下水利用率、供水模数、需水模数、重复利用率、人均供水量、生态环境用水率[1-2]。这一套指标体系已经应用到挠力河流域,并采用主成分分析法、集对分析法进行计算评价[3]。但诸多实践表明,承载力综合评价所得到的结果不尽合理,且地下水承载力综合评价结果的合理性尚未得到有效验证,由此就出现了综合评价指标体系及评价等级合理性的问题[4-5]。

地下水承载力综合评价是否合理,目前主要从评价指标的选取、指标体系的构建、评价方法的选择等方面控制[5],但从地下水承载力综合评价结果合理性角度,目前学者只是对评价结论进行分析,发表的相关研究论文较少。评价结果的合理性,直接影响地下水资源承载能力真实性;因此,本文以水文实测数据为基础,利用较成熟的地下水数值模拟软件MODFLOW,对挠力河流域农灌区地下水进行建模,评价各灌区开采潜力,并以所得结果作为参考,分析地下水承载力综合评价结果的合理性,以期促进本地区相关指标体系的构建与优化。

1 研究区概况

挠力河流域由挠力河不断迁移和泛滥形成,位于黑龙江省东北区域[3],流域总面积约2.5万 km2,南部包括宝清县行政区全部,西、北部包括友谊、富锦、集贤、双鸭山等县市行政区部分区域,东部包括饶河县的部分区域[3];并建有红兴隆和建三江2个农垦管理局,流域包括东泄总灌区、锦西灌区、八五三灌区、大兴灌区、龙头桥灌区、锦南灌区、友谊西部灌区等12个灌区及其他灌区的部分区域[3]。

该流域地处中纬度,为大陆季风气候,雨热同期,年均气温1℃～4℃,年降水量500～600 mm,集中在6-9月。流域内优势河流主要有挠力河和七星河;挠力河全长约596 km,流域总面积24 863 km2,其中挠力河干流流域面积为13 084 km2,支流七星河流域面积为10 815 km2,其他流域面积964 km2;挠力河流域整体为条状,支流基本呈现叶状。该流域因地形地势和水文地质影响,形成了较大面积的沼泽湿地。

挠力河流域整体呈现西南向东北倾的特点,地形相对简单;受构造及河流的影响,根据地貌基本成因类型,本区可划分为堆积地形(I1、I2)、剥蚀堆积(II1)、火山堆积(II2)及侵蚀剥蚀地形(III1、III2)四种。平原上除零星分布孤山和残丘,高度多在500 m以下,主要由古生代、中生代页岩,中酸性火山岩和和花岗岩所构成;多数区域在复合的褶皱基底上堆积有千米以上的中、新生代沉积盖层,第四系厚度一般为120～200 m,最厚可达220 m;大部分孔隙水含水层为砂和砂砾石层,单井出水量一般为150～300 m3/h,适合农业开采。此次模拟该流域第四系孔隙含水层系统,由于第四系较厚,存在4～5个含水层,但各含水层间无明显隔水层,且弱透水层不连续,形成统一大厚度含水层[6]。

2 地下水数值模拟

2.1 水文地质概念模型

2.1.1 水文地质条件概化

研究区域以挠力河流域灌区为核心,向西、向北扩到挠力河与松花江的分水岭,向南、向东分别延伸到平原区与山丘区的交界线;研究区较流域面积小,主要为第四系覆盖区,总面积约为1.40万 km2。

挠力河流域表层广泛覆盖有3～20 m的亚黏土,因此地下水具弱承压性(水头6～7 m)。含水层为第四系松散岩类砂及砂砾石层。挠力河流域水平方向上由山前台地断续分布的第四系微孔隙裂隙水到挠力河流域平原地区的第四系孔隙弱承压水逐步过渡,松散岩类含水层逐渐变厚,粒径逐渐变细。

依据本区域地质情况和实测水文资料,建立适合本研究区的水文地质概念模型;该模型在垂向上将含水层分为2层,即上层为潜水含水层,下层为承压含水层。上层含水层(潜水层)岩性主要是亚黏土,在挠力河流域分布广泛;下层含水层(承压层)主要为砂及砂砾石。含水介质概化成非均质、各向同性。

研究区AB界线为河流流出区,侧向流出研究区,以侧向径流为主;BC和CD界线为挠力河与松花江的流域分水岭;DE和EA界线为山区与平原区的分界线,边缘山地侧向径流补给平原区(图1);因此,整个区域的侧向边界均概化成第三类边界。挠力河、七星河在区内流量较大,与第四系松散岩类含水系统水量交换广且大,可概化为河流边界。研究区上部以潜水面为界,经由该边界浅层地下水与外部环境产生水量交换;底部以第四系松散层底板为界,概化为隔水边界[7]。挠力河流域松散岩类含水层分布广且厚度大,符合达西定律,概化为层流。潜水和承压含水层之间有水量交换,概化为三维非稳定流。研究区的地下水补给项可概化为4项:大气降水入渗、灌溉回渗(井灌和渠灌回灌、渠系回渗)、基岩区地下水径向径流补给和河流渗漏补给[8];地下水排泄项也可以概化为4项:蒸发量、人工开采量(包括工业、农业和生活用水开采量)、河流排泄量、侧向径流排泄量[8]。

图1 研究区水文地质概念模型

2.1.2 水文地质参数分区

主要的水文地质参数为渗透系数与给水度/贮水率。根据挠力河流域地貌类型和水文地质特点,将模拟区潜水含水层渗透系数和给水度划分为2个参数分区,深层含水层渗透系数和贮水率分为8个参数分区(图2)。本次研究主要依据包气带岩性对降水入渗系数以及极限蒸发深度进行分区(图3)。

图2 潜水含水层(a)和承压含水层(b)水文地质参数分区图

图3 包气带岩性分区及入渗、蒸发参数分区图

2.2 地下水流数学模型

(1)

2.3 模型求解

采用Visual Modflow软件系统,对挠力河流域地下水数学模型进行求解,该软件具有强大计算功能以及良好的图形可视界面功能,得到国际同行的广泛认可[9]。

2.3.1 区域剖分和时空离散

根据地下水位动态监测数据的时空分布规律,识别期设为2012年1月1日-2012年12月31日,验证期为2013年1月1日-2013年12月31日,以5 d为1个时段,共69个时段;该时段经历了枯、平、丰水期,能更好的反映含水层系统的特征。

计算区采用矩形网格剖分法,设定网格剖分尺寸为1 385 m×1 385 m,单个计算单元面积为1.9 km2。

2.3.2 含水层参数

根据挠力河流域实际水文地质情况及既往研究成果,按照图3所示含水层分区,设定模型渗透系数与给水度/贮水率,根据各单元含水层岩性特征赋各参数分区的初始值(表1)。按照图4所示降水入渗系数、灌溉回归系数以及极限蒸发深度分区进行赋初值(表2)。

表1 含水层水文地质参数初值

表2 大气降水入渗系数以及极限蒸发深度初始值

图4 模拟值与观测值拟合散点图(a)和模拟误差分布图(b)

2.3.3 模型识别和验证

模型识别是将所建立模型的边界条件、源汇项、水力特征、水文地质参数等与实际情况相吻合,通过不断调整模型参数,各时段的地下水位与实测的地下水位拟合度较好(图4a),两者之间的相关系数可达0.99;误差分布直方图基本呈现正态分布,拟合误差最小点出现的频率位于正态分布曲线的中心部位(图4b)。

经识别后的含水层渗透系数、给水度的参数分区值如表3;大气降水入渗系数分区、大气降水入渗系数分区值如图5、表4所示;极限蒸发深度分区及赋值与初始一致(表2)。

表3 含水层渗透系数和给水度初值

表4 大气降水入渗系数识别值

图5 大气降水入渗系数识别分区图

将识别所得的水文地质参数输入到模型中,验证期为2013年1月1日-2013年12月31日,以5天为1个时段,检验下水位计算值与观测值的拟合情况,结果显示地下水位拟合效果良好,计算值与实测值的相关系数达到0.993(图6a);误差分布直方图呈现正态分布,拟合误差最小点出现的频率位于正态分布曲线的中心部位(图6b)。

图6 模拟值与观测值拟合散点图(a)和模拟误差分布图(b)

经过模型的识别与验证,所建模型的边界条件、源汇项、水力特征、水文地质参数等与实际情况拟合效果良好,能够正确反映区域水文地质条件以及地下水运动规律,可以采用此模型对研究区地下水动态规律开展模拟研究。

3 地下水资源开采潜力评价

Modflow软件中的Zone Budget模块可对整个区域进行水量均衡分析,也可对局部区域进行水量均衡分析;利用Zone Budget模块可统计出模型中的区域入渗补给量、边界补给量、河流补给量以及排泄量等数据。利用Zone Budget模块,将研究区内的各个灌区设为具体的水量均衡区进行统计分析(图7)。

图7 水量均衡统计分区示意图

模型中地下水总补给量包括入渗补给量、河道渗漏补给量、边界侧向径流补给量、区间侧向径流补给量等。即:

Q总补= Q渗补+ Q河补+ Q边补+ Q侧补

(2)

式中:Q总补为地下水总补给量,104m3/a;Q渗补为包括:降水入渗补给量、井灌回归补给量、渠系渗漏补给量、渠灌田间入渗补给量,104m3/a;Q河补为河道渗漏补给量,104m3/a;Q边补为由通用水头边界(GHB)计算出的边界侧向径流补给量,104m3/a;Q侧补为相邻灌区之间侧向径流补给量,104m3/a。

灌区现状年地下水可开采量的计算采用可开采系数法,即:

Q可采= ρ· Q总补

(3)

式中:Q可采为地下水可开采量,104m3/a;Q总补为地下水总补给量,104m3/a。ρ为可开采系数(ρ≤1),由于本流域地下水富水性较好,因此对于单井单位降深出水量>20 m3/h·m,且地下水埋深大、水位连年下降的超采区,ρ值的取值范围为0.8～1.0;对单井单位降深出水量在5～10 m3/h·m,地下水埋深大、实际开采程度较高地区或地下水埋深较小、实际开采程度较低地区,ρ值的取值范围为0.6～0.8[10-12]。

经计算可得挠力河流域灌区现状年地下水总补给量、地下水可开采量、资源量及可开采量,详细计算见表5。

表5 灌区现状年地下水总补给量、可开采量及资源量计算表

基于表5的结果,参考杨国强(2012)在地下水超采指标体系建立中,利用开采系数法对研究区地下水开发利用情况划分等级,即:小于95%为潜力区,95%～105%为补给平衡区,105%～120%为一般超采区,大于120%为严重超采区[13]。对挠力河地下水开发利用情况划分为5个等级,即:Ⅰ级开采率低于50%(潜力较大);Ⅱ级50%～80%(潜力中等);Ⅲ级80%～100%(潜力较小);Ⅳ级100%～120%(一般度超采);Ⅴ级大于120%(重度超采)。

由表6可知,挠力河流域灌区现状年总的可开采量为1.44×109m3/a,总的开采量为1.63×109m3/a,地下水开采率为113.1%,灌区总体呈超采状态。但各灌区地下水资源的承载力不同:(1)蛤蟆通灌区、七里沁灌区、尖山子灌区、三环泡灌区、龙头桥灌区开采率低于50%,属开采潜力较大的灌区,尚具有较大的农业灌溉潜力;(2)集贤小黄河灌区、幸福灌区开采率在50%～80%之间,属开采潜力中等,还可支持一定规模的农业灌溉面积发展;(3)富锦灌区地下水开采率为83.5%,开采率在80%～100%之间,具有较小的开采潜力,如果继续大规模开采,将有超采的可能。(4)二九一南部灌区、友谊西部灌区开采率在100%～120%之间,属一般超采区。灌区地下水开采量已经超过该灌区地下水可开采量,如果灌区地下水利用量不加调节则地下水位将持续下降。(5)锦西灌区、五九七灌区、东泄总灌区、锦南灌区、大兴灌区、八五三灌区地下水开采率大于120%,属重度超采,部分灌区已出现地下水降落漏斗,若不加调节将可能引发严重的环境地质问题。

表6 灌区地下水开发利用程度分析

4 承载力综合评价结果合理性分析

地下水资源承载力强弱程度是区域内地下水资源系统与社会、环境综合作用的结果,直观地表现在地下水开采潜力上,区域地下水承载力强则表现为地下水开采潜力大,反之则表现为地下水开采潜力小;根据这一特性,本文采用已有的挠力河流域灌区地下水承载力综合评价结果[3],其结果见表8;并以模拟所得出的结果作为参考,评价其合理性。

因两个体系的分级不同,因此为了便于两者的比较,根据两者等级划分标准,给出了两者的对比标准(表7),当两者一致度达到80%以上时即可接受,当低于70%时需要重新修订指标体系[14]。参考已有的地下水承载力综合评价结果,将地下水承载力综合评价结果和地下水开采潜力分析结果汇总于表8。

表7 一致性判别表

表8 地下水承载力综合评价结果合理性分析

由表8可知,从排序上看,除二九一南部灌区、富锦灌区、五九七灌区和三环泡灌区外,两种方法所得的评价结果都比较相近;但从评价等级上看两者相差较大,二九一南部灌区、锦西灌区、锦南灌区、友谊西部灌区、五九七灌区地下水承载力为中或高等水平,但运用数值模拟分析的实际情况为严重超采区,开采潜力很弱。此外,地下水承载力综合评价结果在等级划分上不够细化,如:集贤小黄河灌区、幸福灌区地下水承载力为高等,但该灌区实际开采潜力为中等水平,这就导致评价结果出现偏差,影响评价结果的合理性。从一致性判断上分析,两者的一致性为68.7%,一致性低,需要对指标体系重新修订。

分析其中的原因主要包括4个方面:(1)挠力河流域乃至整个三江平原区域内耕地面积多,人口少,以人口量作为区域衡量指标不合理,所选人均供水量指标欠妥。(2)区域地下水资源可利用量计算较复杂,此外,如果已知区域地下水开发利用程度,也就已知区域地下水整体开采情况,再以地下水水量为基础选取指标建立指标体系,进而评价地下水资源承载力是不合理的。(3)挠力河流域灌区内第二、三产业几乎不存在,而农业较发达,在后续选取指标时,应以“地下水—农业”二者之间关系为主。(4)关于地下水资源量,其计算也比较复杂,对于挠力河流域各灌区灌溉期集中,且灌区地下水主要补给来源为降水入渗补给,可以降水入渗量作为地下水资源量。

综合以上分析结果来看本次地下水资源承载力评价所采用的指标体系尚存在一定的不合理性。因此,指标选取需要仔细甄别,指标分级标准也需要细化,以便构建一套合理的地下水资源承载力评价指标体系。