泰安旧县水源地地下水位动态特征及可开采量研究.pdf

泰安旧县水源地地下水位动态特征及可开采量研究魏凯1,2，王延岭1,2，赵志伟1,2，吴亚楠1,2，翟代廷1,2，闫佰忠3（1.山东省地质矿产勘查开发局第五地质大队,山东泰安271000；2.山东省地矿局岩溶塌陷防治重点实验室/泰安市岩溶塌陷防治重点实验室,山东泰安271000；3.河北地质大学水资源与环境学院,河北石家庄050031）摘要：泰安市城区水资源供需矛盾突出，地下水资源开采引发了岩溶塌陷等地质环境问题，亟需查明岩溶塌陷生态水位约束下水源地的允许开采量，确定水源地最优取水方案。以旧县水源地19802021 年的地下水位、年降水量、地下水开采量、岩溶塌陷等数据为基础，探讨人类活动对水源地地下水位的演变的影响因素，分析发生岩溶塌陷的地下水临界水位。基于 Modflow-GWM 软件构建泰安城区旧县岩溶水系统流动与管理耦合模拟模型，探讨防止岩溶塌陷发生的地下水允许开采的资源量。结果表明：（1）孔隙水水位年内动态受降水影响明显，呈现“降水补给”型特点；岩溶水水位变化呈“枯低丰高”的特征，水位动态变化属降水入渗开采型；（2）旧县水源地地下水位自19801990 年呈现大幅下降，19902003 年岩溶地下水位基本呈现波动下降，2004 年后水源地地下水位上升较为明显；（3）研究水位动态与岩溶塌陷关系得出，防止发生岩溶塌陷的水源地临界水位为 108m，处于岩溶含水层顶板以上 2m；（4）通过地下水管理模型，确定在临界水位时模拟区岩溶水可开采量为 8.28.5万 m3d1，其中旧县水源地可开采资源量为 3.23.5 万 m3d1。关键词：地下水水位；人类活动影响；人工开采；岩溶塌陷；地下水优化模型中图分类号：P641.8文献标识码：A文章编号:10014810（2023）05094016开放科学（资源服务）标识码（OSID）：0引言在地表水资源相对匮乏的我国北方地区，地下水是维持人们正常生产生活的重要水源1。但是，由于地下水的大规模集中开采，导致了地下水降落漏斗、地面沉降和岩溶塌陷等一系列生态环境问题。目前，我国北方地区针对地下水超采问题采取了一些治理措施，地下水超采问题得到了明显改善2。但是局部地区地下水的严格禁采也引发了诸如地下水浸没、沼泽化、盐碱化、资源利用率低等诸多问题3。如何在保障地下水生态环境健康的前提下，最优化开发利用地下水资源，成为当前我国北方地区生态文明建设的重大问题。泰安市地处鲁中山区，工农业生产发达，水资源相对匮乏，而城区水资源供需矛盾尤为突出，特别是极端干旱年份，城市供水安全受到严重威胁。泰安市城区供水主要依靠黄前水库地表水源以及旧县、大汶口两处地下水源供给，地下水源供水占比达40%以上。由于泰安市城市规模的不断扩大，城区水资源供求矛盾日趋突出，地表水源供水能力有限，基金项目：泰安市科技创新发展项目（2020NS291）；国家自然科学基金项目（42002251）第一作者简介：魏凯（1986），男，高级工程师，主要从事水工环地质工作。E-mail：。通信作者：王延岭（1978），男，高级工程师，主要从事水工环地质工作。E-mail：。收稿日期：20230420第42卷第5期中国岩溶Vol.42No.52023年10月CARSOLOGICASINICAOct.2023魏凯，王延岭，赵志伟，等.泰安旧县水源地地下水位动态特征及可开采量研究J.中国岩溶，2023，42（5）：940-955.DOI：10.11932/karst20230507地下水资源开采量逐年增大，引发了岩溶塌陷等地质环境问题。在泰安市城区的中南部地下水开采强度较大的地区形成了以城区水源地塌陷区、旧县塌陷区和羊娄塌陷区为代表的岩溶塌陷带，面积达到200km24。为了防止岩溶塌陷的进一步加剧，泰安市对城区水源地和旧县水源地分别进行了限采，其中城区水源地供水井 1993 年后全部关停；旧县水源地 2003 年后逐步进行限采5。地下水水源地的限采对泰安市的供水安全产生了一定的影响，因此，如何确定岩溶塌陷生态水位约束下水源地的允许开采量，进而提出旧县水源地可开采资源量和最优取水方案，保障水源地地下水资源合理开发和管理规划至关重要。目前，对于地下水水源地可开采量的计算主要采用水均衡法和数值模拟法。周博文等6从水文地质条件入手，采用水均衡法对额仁淖尔水源地靶区地下水可开采潜力进行了分析评价。李若怡等7考虑植被、湖泊的生态保护目标约束，通过地下水流模型重新评价了水源地的可开采资源量，优化了开采井布局方案。但是这两种方法只能解决在给定方案前提下的预测问题，不能确定哪个方案能获得最优结果；而通过数值模型和管理模型耦合进行优化，可给出满足设定目标函数和约束条件下的最优开采方案8。对于水源地管理优化模型，合理约束地下水水位的确定对于优化模型是至关重要的。Alfaro 等9以地下水天然水位作为约束探讨了南苏丹河谷区的地下水超采问题。马雄德等10探讨了面向生态的矿区地下水位控制阈值。付晓刚等11在探讨羊庄盆地地下水可开采量时，将盆地未发生岩溶塌陷的地下水位历史最低值作为约束水位。可见将地下水流数值模型与管理模型耦合，确定地下水约束水位是分析水源地合理开采量的关键。本文以旧县水源地 19802021 年的地下水位、年降水量、地下水开采量、岩溶塌陷等数据为基础，探讨了人类活动影响下泰安市旧县水源地地下水位的演变特征及影响因素；并结合岩溶塌陷和地下水位的关系，确定合理的防治岩溶塌陷的临界水位，基于 Modflow-GWM 软件构建了水源地的地下水流动与管理耦合模拟模型，对防止岩溶塌陷的临界地下水位及允许开采资源量进行探究，确定旧县水源地的最优取水方案，以期能为岩溶地下水水源地合理开发和管理规划提供科学依据。1研究区概况旧县水源地位于泰安城区旧县岩溶水系统的东南部，岩溶水系统北至泰安城区北部环山路，南至牟汶河以南，西至泰安城区西 G104 国道，东至邱家店镇埠阳庄，东西长约 25km，南北宽约 12km，地理 坐 标 为：东经 11703 43 11718 16，北 纬360519361222，研究区属温带大陆性半湿润季风气候区，四季分明，春旱多风，夏热多雨，秋旱少雨，冬寒少雪，季节性干旱严重。根据多年（19502020 年）统计资料，气温、降水量和蒸发量情况如下：平均气温 1113，多年平均降水量 765.3mm，年最大 降 水量 1 571.7mm（1964 年），年 最 小 降 水 量263.1mm（1989 年），降雨量时空变化较大，主要集中在汛期 79 月份，约占全年降水量的 66%；年总蒸发量 1664.21927.0mm。研究区为汶河水系，区内地表水系十分发育，以季节性河流为主，其主干河流为牟汶河，由北东向南西穿过研究区东部。研究区在区域上位于泰莱盆地西南隅，北部为泰山山脉，南部为徂徕山脉，两山脉均为山势陡峻、切割强烈的中低山。整个地势由东向西倾斜，南、北两面又向中部盆地倾斜。地貌类型为山间冲洪积平原，地势总体东北高，西南低，地面标高一般 120140m，坡度小于 5。区域地层发育有太古界前震旦系泰山群变质岩，下古生界寒武系、奥陶系碳酸盐岩、页岩，新生界古近系砾岩、泥岩、泥灰岩、黏土岩及第四纪松散岩等（图 1）。泰安市城区旧县岩溶水含水系统内的含水岩组包括第四系松散岩类孔隙含水岩组和寒武奥陶系碳酸盐岩类裂隙岩溶含水岩组两种类型。第四系孔隙含水层为牟汶河沿岸冲积及山前倾斜平原冲洪积层，岩性为中、粗砂夹卵砾石，一般分选性较好，除河床内直接裸露外，皆上覆有一定厚度的砂质黏土，此盖层具有一定的渗透能力，为降水垂直渗入造成了有利条件。第四系含水层与隐伏寒武系和奥陶系地层之间无明显隔水层，彼此间存在水力联系。第四系孔隙含水层主要接受地表水以及大气降雨的补给，径流方向自北向南指向牟汶河，与地形坡度保持一致，排泄方式为补给地表水、蒸发、越流补给以及第42卷第5期魏凯等：泰安旧县水源地地下水位动态特征及可开采量研究941人工开采。第四系松散岩类含水岩组包括冲洪积层及冲积层两种类型，岩性的差异造成富水性不均。其中第四系冲积孔隙水含水层呈条带状沿牟汶河分布，含水层上部为浅黄、褐黄色中粗砂，下部为粗砂夹卵砾石，磨圆度较好，单井涌水量在 10003000m3d1之间。第四系冲洪积孔隙含水层主要分布于牟汶河以北山前倾斜平原，含水层上部为细砂，中部为粗砂，下部为砂砾石，磨圆度较差，单井涌水量小于 3000m3d1。碳酸盐岩裂隙岩溶含水岩组是本次重点关注的含水岩组。含水层主要为寒武系上统凤山组及奥陶系下统和中统二、三段之石灰岩，由南向北大致呈北西走向分布，绝大部分被第四系堆积物覆盖。灰岩中裂隙岩溶较发育。基岩表面裂隙岩溶的发育有利于岩溶水的补给，深部裂隙岩溶的发育有利于岩溶水的运动和富集。区内岩溶水主要接受孔隙水越流补给、牟汶河渗漏补给，主要排泄方式主要是人工开采。其中旧县水源地取水层位主要为碳酸盐岩裂隙岩溶含水岩组，旧县水源地自 2003 年后逐步进行限采。目前，旧县水源地集中量 1.8 万 m3d1左右，另外区内农村分散式生活供水和农用灌溉井有近350 口，日均开采量 5 万 m31213。2地下水位动态及岩溶塌陷关系2.1地下水位动态特征孔隙水年内动态受降水影响明显，年最低水位一般出现在 6 月下旬，年最高水位一般出现在雨季的 8 月末至 9 月末，呈现“降水补给”型特点。水位年变幅维持 15m。旧县监测点受降水与旧县水源地地下水开采双重影响，水位波动较径流区稍大，19902002 年间，孔隙水水位小幅下降，20032005年受降水量大影响，水位快速回升，2006 开始，水位又 呈 下 降 趋 势，至 2019 年 出 现 历 史 最 低 水 位108.47m（2019 年 6 月），受年度降水量大影响，2020年地下水水位有所上升，年末水位 130.72m。该点多年水位动态整体呈下降趋势，水位年际变幅范围在 0.724.49m（图 2，图 3）。旧县水源地岩溶地下水年内动态受第四系孔隙水排泄受阻影响，年内最高水位一般略微滞后于集中 910 月，最低水位则出现在 56 月，年内地下水位变幅 5m 左右。水位变化呈“枯低丰高”的特征，水位动态变化属降水入渗开采型。自 1990 年以来，旧县水源地负担泰安供水，加上农业灌溉影响，岩溶2 km 024 km沂河组临沂组黑土湖组大站组朱家沟组三山子组炒米店组崮山组张夏组馒头组朱砂洞组地质界线断层研究区范围马家沟群图1旧县水源地区域地质图Fig.1GeologyofthegroundwatersourceareainJiuxiancounty942中国岩溶2023年地下水水位变化受开采量影响较大，水位年际动态呈现出 2 个阶段特征。19902003 年岩溶地下水位基本呈现波动下降；2003 年限采后岩溶地下水位抬升后基本达到稳定，从 2003 年的 111.20m 升至2021 年水位 113.98m（图 4，图 5）。2.2岩溶塌陷与地下水位关系根据监测资料以及前人研究14。1982 年水源地开始向泰安城区供水，初期供水量 1.382.73 万m3d1，水位由开采前的 114m（埋深 46m）下降至1983 年的 108m（埋深 12m）左右，开始出现塌陷；此后水位有所回升，但随着开采量的增大，水位整体呈现下降趋势，特别是自 1987 年开始，供水量增加至4.63 万 m3d1，水位急剧下降，由 108m（埋深约 12m）左右迅速下降至 1990 年的 96m（埋深 24m）左右，自 1988 年开始，旧县水源地进入岩溶塌陷高发期，此后至 1994 年，由于开采量有所减少（4 万 m3d1），水位有所回升，在 96108m（埋深 1224m）之间波动，处于基岩面附近，这期间也是区内岩溶塌陷的高发期，塌陷区范围也由旧县向北扩展到苑庄羊娄一带；1995 年之后，区内岩溶塌陷的强度有所降低，但仍时有发生，至 2001 年为区内塌陷的发展期，期间水位总体呈下降趋势，并在基岩面上下波动，水位标高 95110m（埋深 1025m）；2002 年为大旱之年，全年降水量不足 300mm，水源地自来水厂开采量达到历史最高的 5 万 m3d1，同时农灌开采量激增，年最低水位下降至 91m（埋深 29m）左右，2003 年上半年持续高强度的开采，水位在 8698m（埋深 2234m）之间波动，高强度的开采和持续的低水位在 2013 年降水量降雨量/mm地下水位标高/m孔隙地下水水位时间/月2021/012021/022021/032021/042021/052021/062021/072021/082021/092021/102021/112021/12400350300250200100500150117116115114113112111110109108图22021 年旧县孔隙水 296 监测孔年内地下水位动态降雨量曲线Fig.2Annualgroundwaterleveldynamics-precipitationcurveofNo.296monitoringholeforobservingporewaterinJiuxiancountyin2021降水量孔隙地下水水位时间/月1990/011991/091993/051995/011996/091998/052000/012001/092003/052005/012006/092008/052010/012011/092013/052015/012016/092018/052020/012021/096005004002001000300降雨量/mm孔隙地下水位标高/m122120118116114112110108106104102图3旧县孔隙水 296 监测孔 19902021 年年际地下水位动态降雨量曲线Fig.3Interannualgroundwaterleveldynamics-precipitationcurveofNo.296monitoringholeforobservingporewaterinJiuxiancountyfrom1990to2021降水量降雨量/mm地下水位标高/m岩溶地下水水位时间/月2021/012021/022021/032021/042021/052021/062021/072021/082021/092021/102021/112021/12400350300250200100500150117116115114113112111110109108图42021 年旧县岩溶水 225 监测孔年内地下水位动态降雨量曲线Fig.4Annualgroundwaterleveldynamics-precipitationcurveofNo.256monitoringholeforobservingkarstwaterinJiuxiancountyin2021降水量岩溶地下水水位时间/月1991/011992/031993/051994/071995/091996/111998/011999/032000/052001/072002/092003/112005/012006/032007/052008/072009/092010/112012/012013/032014/052015/072016/092017/112019/012020/032021/05降雨量/mm地下水位标高/m50045040035030025020015010050013012011010090807060图5旧县水源地监测点岩溶水 19912021 年年际地下水位动态降雨量曲线Fig.5Interannualgroundwaterleveldynamics-precipitationcurveofkarstwaterattheobservationpointofgroundwatersourceareainJiuxiancountyfrom1990to2021第42卷第5期魏凯等：泰安旧县水源地地下水位动态特征及可开采量研究943雨季引发了大范围、高强度的岩溶塌陷地质灾害，2003 年下半年水源地开采量减少至 3.8 万 m3d1，农灌井开采量急剧减少，水位迅速回升至 115m（埋深5m）左右，2004 年水位在 105116m 之间波动，又发生多次塌陷，但规模较小；2005 年至今，水位总体呈现上升趋势，在 106116m 之间波动，位于基岩面以上，期间 2006 年发生塌陷一处，2013 年和 2014 年旧县村发生房屋多处房屋斑裂，未见形成新的塌陷坑（图 6）。据此，确定旧县水源地地下水位生态约束水位为 108m。地层年度代号岩石名称 柱状图坐标轴分类耕植土黏土中粗砂房屋斑裂塌陷数量开采量降水量塌陷数量/处石灰岩Q8590952 000051015202466543281012141618130120100806040降水量/cm开采量/万 m3d1H/D(m)图6旧县水源地岩溶塌陷发育特征与水位关系（吴亚楠，2017 年）Fig.6RelationshipbetweenthedevelopmentcharacteristicsofkarstcollapseandgroundwaterlevelinthegroundwatersourceareaofJiuxiancounty(WUYanan,2017)3地下水流动模型构建3.1水文地质概念模型3.1.1模拟区范围区内孔隙水模拟范围基本与岩溶水模拟范围相同，但是在西南部南刑家寨朱家埠一带，寒武系地层直接出露地表；此外，在模拟区的桂林官庄西颜张一带和模拟区的东部徂徕镇鲁家庄禄官路一带寒武系地层也同样出露地表，因此孔隙水模拟范围内不包括这些地区，孔隙水模拟范围面积为 130.986km2。区内岩溶水模拟范围北部边界以泰山断裂为界，断裂北侧为古老变质岩，富水性和透水性较弱，可概化为弱透水边界。东北部边界：中北段以岱道庵断裂为界，东侧为厚度大、透水性弱的古近纪黏土岩分布，西侧为奥陶系灰岩，为相对隔水边界；中南段边界以埠阳断裂为界，该断裂为张扭性断裂，属透水弱透水边界。西部、南部边界为寒武纪岩层与古老变质岩接触带，透水性弱，属弱透水边界。由于东部的边家庄牛角一带缺少岩溶水统测水位和岩溶水观测孔，因此东部边界以逯家庄鲁家庄为界，整个岩溶水模拟范围面积为 155.875km2（图 7）。3.1.2含水层结构概化本次模拟含水层概化为两层结构，第一层为孔隙潜水含水层，第二层为岩溶含水层。第一层孔隙含水层接受大气降水入渗补给，第一层与第二层可通过越流发生水力联系。区内地下水补给来源主要为大气降水、外围地下水的侧向补给、南部牟汶河的渗漏补给、农灌回渗补给等，两层水之间存在越流补给，尤944中国岩溶2023年其是在岩溶天窗部位。排泄途径为人工开采，尤其是岩溶水，区内包含岩溶水水源地，用于城镇生活用水、工业用水，另外还有居民自备井百余口，用于农灌或生活用水。丰水期地下水位高于地表水时，地下水也会排泄于地表水。岩溶含水层厚度约 120m，底部边界概化为不透水边界，与下伏含水层无水力联系。3.1.3边界条件概化垂向边界：含水层概化为两层结构，第一层孔隙潜水含水层，第二层为岩溶含水层。第一层孔隙含水层接受大气降水入渗补给，第一层与第二层可通过越流发生水力联系。区内地下水补给来源主要为大气降水、外围地下水的侧向补给、南部牟汶河的渗漏补给、农灌回渗补给等，两层水之间存在越流补给，尤其是在岩溶天窗部位。排泄途径主要为人工开采，尤其是岩溶水，区内包含三处岩溶水水源地，用于城镇生活用水、工业用水，另外还有居民自备井百余眼，用于农灌或生活用水。丰水期地下水位高于地表水时，地下水也会排泄于地表水。岩溶含水层以埋深约 800m作为底边界，为不透水边界，与下伏含水层无水力联系。平面边界：孔隙水整体流向为由西北向东南流动，模型中将西北和东部边界概化为第二类流量边界，其他边界概化为零通量边界。岩溶含水层由于水源地的开采，整体流动向旧县水源地漏斗区汇集，北部边界以泰山断裂为界，概化为弱透水边界。东北部边界：中北段以岱道庵断裂为界，为相对隔水边界，概化为零通量边界；中南段边界以埠阳断裂为界，属透水弱透水边界，概化为弱透水边界。西部、南部边界为寒武纪岩层与古老变质岩接触带，透水性弱，概化为弱透水边界。3.1.4源汇项处理模拟区内孔隙含水层主要接受大气降雨入渗补给、地下水侧向径流补给和牟汶河渗漏补给，排泄方式主要是越流排泄和人工开采；岩溶含水层主要接受越流补给、地下水侧向径流补给和牟汶河渗漏补给，排泄方式主要是人工开采。其中大气降雨入渗补给条件的不均匀性在模型中用降雨入渗系数分区概化处理；孔隙水和岩溶水农业灌溉开采利用面状开采处理；水源地集中开采利用 wells 模块处理；牟汶河渗漏补给利用 river 模块处理。3.2数学模型及离散化研究区地下水补排量和水位动态是随时间变化的，表现出非稳定流特性；孔隙含水岩组的岩性在空城市乡镇岩溶水观测孔隔水边界弱透水边界模拟区范围断裂带河流024 km图7模拟范围及边界条件概化Fig.7Generalizationofsimulationscopeandboundaryconditions第42卷第5期魏凯等：泰安旧县水源地地下水位动态特征及可开采量研究945间上各地存在一定差异，反映了系统的非均质性，垂向由于冲积的成层性和后期的压缩作用使之与水平方向异性；岩溶含水层组也具有非均质性，同时由于受构造影响，岩溶发育体现出一定的各向异性特征。因此，总的来说本地区地下水系统可当非均质、各向异性、三维非稳定流模型来进行研究。SHx=x(KxHx)+y(KyHy)+z(KzHz)+x,y,z,t 0Ht=Kx(Hx)2+Ky(Hy)2+Kz(Hz)2Hz(Kz+p)+px,y,z 2,t 0h(x,y,z,t)|t=0=H0 x,y,z,t 0Knhrn?2=q(x,y,t)x,y,z 2,t 0式中：渗流区域；H含水层水位标高（m）；H0含水层初始水位（m）；Kx、Ky、Kz-分别为 x、y、z 方向的渗透系数（md1）；S承压含水层储水率；孔隙含水层给水度；含水层的源汇项（d1）；0渗流区域的上边界，即地下水的自由表面；2渗流区域的第二类边界，包括含水层隔水底边界和渗流区域的侧向流量边界；p孔隙水的蒸发和降水补给等（d1）；nr边界面的法线方向；Kn边界面法向方向的渗透系数（md1）；q（x，y，z，t）定义为二类边界的单位面积流量（m3d1）。空间离散：将模拟区剖分为 300m400m 的规则矩形网格，共剖分单元格数量 240000 个，其中有效单元格数 108999 个。3.3模型识别验证根据数值模拟时期的确定原则，将 2016 年 6 月作为水位模拟的初始时刻（图 8，图 9）。利用 2016年 6 月2021 年 12 月的水位、源汇项数据资料进行模拟，对模型进行识别验证，其中 2016 年 6 月2019 年 12 月为模型识别期，2020 年 1 月2021 年12 月为模型验证期。应力期为一个月，每个应力期分为三个时间步长，共 67 个时间步长。孔隙水水位动态观测孔共 4 个，岩溶水水位动态观测孔共8 个。024 km孔隙水地下水位模拟区范围城市乡镇河流断裂带图82016 年 6 月孔隙水等水位线图Fig.8GroundwaterlevelcontourmapofporewaterinJune,2016946中国岩溶2023年3.4模型识别验证识别验证综合利用孔隙水和岩溶水地下水位动态长观孔以及模拟区流场进行拟合。经过对模型输出结果分析，识别验证期的地下水位动态观测孔水位过程线均拟合较好（图 10，图 11）。在识别验证期孔隙水的平均误差不超过 0.5m，相对误差为 0.843%；岩溶含水层的平均误差不超过 1.0m，相对误差为1.546%（图 12、图 13）。3.5模型识别验证后参数根据模型最终调试结果，识别验证后降水入渗系数分区见表 1，孔隙水渗透系数和给水度分区见表 2，岩溶水渗透系数及储水系数分区见表 3。其中降水入渗系数共分 8 个区，值区间为 0.120.30；孔隙水渗透系数共分为 13 个区；岩溶水渗透系数共分为9 个区（图 14 至图 17）。4地下水优化管理模型构建4.1最优化管理模型在地下水流动数值模型的基础上（图 14），利用GroundwaterManagement（GWM）构建地下水管理模024 km岩溶水地下水位模拟区范围城市乡镇河流断裂带图92016 年 6 月岩溶水等水位线图Fig.9GroundwaterlevelcontourmapofkarstwaterinJune,2016201712512020182019202020212022地下水位/m272 孔隙水观测孔Trans,Head272272-Observed时间/a图10识别验证期孔隙水 QS8 观测孔实测地下水位与计算地下水位拟合曲线Fig.10FittingcurveofmeasuredandcalculatedgroundwaterlevelofporewateratQS8observationholeinidentificationandverificationperiod201711211411611810811020182019202020212022地下水位/m225 岩溶水观测孔Trans,Head225225-Observed时间/a图11识别验证期岩溶水 225 观测孔实测地下水位与计算地下水位拟合曲线Fig.11FittingcurveofmeasuredandcalculatedgroundwaterlevelsofkarstwateratNo.225observationholeinidentificationandverificationperiod第42卷第5期魏凯等：泰安旧县水源地地下水位动态特征及可开采量研究947型，结合约束条件，在不发生岩溶塌陷的前提下，确定旧县水源地的最大可开采量。因此，本次将旧县水源地 6 眼地下水开采井总开采量最大作为目标函数，把单个水井开采量作为决策变量。考虑到旧县水源地 6 眼水井开采量基本是均分状态，故将各水井开采量均作为不随时间变化的常量进行优化计算，024 km预测地下水位计算地下水位模拟区范围城市乡镇河流断裂带图12识别验证期（2021 年 12 月）孔隙水流场拟合图Fig.12FittingdiagramofporegroundwaterflowfieldduringtheidentificationandvalidationperiodinDecember,2021024 km预测地下水位计算地下水位模拟区范围城市乡镇河流断裂带图13识别验证期（2021 年 12 月）岩溶水流场拟合图Fig.13FittingdiagramofkarstwaterflowfieldduringtheidentificationandvalidationperiodinDecember,2021948中国岩溶2023年据此建立优化模型。目标函数为旧县水源地开采量最大，即：Q总=max6i=1Qi式中：Q总为 6 眼开采井的总开采量；Qi为第 i 眼开采井的开采量。约束条件，水位约束：旧县水源地 6 口开采井在不同时刻的水位均不应低于该处最低约束水位，根据地下水水位动态与岩溶塌陷规律确定；非负约束：决策变量应均为正值，即各开采井开采量均为正值。采用响应矩阵法将前面校正和验证后的地下水流模型与该优化管理模型（包括目标函数、约束条件）耦合起来，通过模型联合运算。4.2优化模型条件的确定4.2.1优化模型中源汇项确定区内地下水的补给主要是大气降雨入渗补给和牟汶河河水渗流补给，因此，未来降雨量的选择对于岩溶水水位预测影响较大，本次采用 19512022 年72 年的降水系列作 P-III 型频率分析，确定丰、平、枯水年降雨量，结合降雨序列的周期性，选择涵盖丰、平、枯周期变化的 15a 降雨序列作为未来 15a 的降雨可能值，平均降水量为 718.2mma1，为一相对偏枯降水序列。由于其他地下水源汇项历年变化不大，仍保持现状。此外，区内分散式生活用水和农业灌溉主要开采岩溶水，这部分开采量仍按照 2020 年保持现状。4.2.2开采约束条件为实现旧县水源地地下水环境持续好转的总目标，防止岩溶塌陷发生，结合 2.2 节地下水水位动态与岩溶塌陷关系，确定旧县水源地岩溶水约束地下水位为 108m。X-axis/103510515520525图14模拟区数值模拟结构示意图Fig.14Numericalsimulationstructureinthesimulatedarea024 km降雨入渗系数分区12345678图15模拟区大气降水入渗系数分区Fig.15Zoningofatmosphericprecipitationinfiltrationcoefficientinthesimulatedarea第42卷第5期魏凯等：泰安旧县水源地地下水位动态特征及可开采量研究949024 km渗透系数和给水度分区12345678910111213图16模拟区孔隙水渗透系数和给水度分区Fig.16Partitionofporegroundwaterpermeabilitycoefficientandwateryieldinthesimulatedarea024 km渗透参数和贮水系数分区1023456789图17模拟区岩溶水渗透系数和储水系数分区Fig.17Partitionofpermeabilitycoefficientandstoragecoefficientofkarstwaterinthesimulatedarea950中国岩溶2023年5结果与讨论使用 GMS 中 的 MODFLOW2000 模 块 建 立MODFLOW 模型，经过识别与验证之后，从中提取GWM 所需的 MODFLOW 数据文件，然后在 GWM软件中输入优化计算所需要的决策变量、目标函数、约束条件以及优化方法，最后利用 GWM 所提供的序列线性规划法来求解该优化模型（图 18 至图 22），从而得出在防止岩溶塌陷前提下，保持模拟区农业灌溉和分散生活岩溶水开采量为 5 万 m3d1时，旧县水源地在集中开采方式下，可开采量为 11681278万 m3年1（即 3.23.5 万 m3d1）（表 4）。6结 论（1）旧县水源地孔隙水水位年内动态受降水影响明显，呈现“降水补给”型特点，多年水位动态整表 1识别验证后模拟区降水入渗系数取值一览表Table1Valuesforprecipitationinfiltrationcoefficientinthesimulatedareaafteridentificationandverification参数分区参数分区参数分区描述10.12山前基岩裸露区20.20丘陵山前地带，岩性为砂质黏土夹砾石粉砂，厚度薄30.30牟汶河河道两侧，岩性为砂含砾石40.14基岩裸露区50.25位于牟汶河西侧，土地类型主要为农耕地、林地，植被覆盖率较高60.22位于泰安城区及以西，岩性主要为粉质黏土夹中粗砂70.22分布于城区水源地周边80.21丘陵山前地带表 2识别验证后模拟区孔隙水渗透系数和给水度赋值Table2Poregroundwaterpermeabilitycoefficientandyieldvalueinthesimulatedareaafteridentificationandverification分区水平K/md1垂向/md1给水度分区水平K/md1垂向K/md1给水度10.150.000100.04083.000.000300.2320.150.000100.04090.850.000150.3032.800.000180.070100.750.000150.3043.000.000300.230113.000.000500.35515.501.500000.200120.500.000100.2560.150.000150.045135.000.000500.3573.000.000300.230表 3识别验证后模拟区岩溶水渗透系数和贮水系数赋值Table3Valuesofpermeabilitycoefficientandstoragecoefficientofkarstwaterinthesimulatedareaafteridentificationandverification分区水平K/md1垂向K/md1弹性储水率/m1015.0104.006.010415.00.505.010424.50.453.010432.50.252.510442.80.283.010552.00.205.510462.30.213.510574.00.405.010584.50.809.510594.00.205.0105第42卷第5期魏凯等：泰安旧县水源地地下水位动态特征及可开采量研究951体呈下降趋势，水位年际变幅范围在 0.724.49m。19902002 年间，孔隙水水位小幅下降，20032005年受降水量大影响，水位快速回升，2006 年开始，水位又呈下降趋势，至 2019 年出现历史最低水位108.47m，2020 年后地下水水位有所上升。（2）旧县水源地岩溶水水位变化呈“枯低丰高”表 4模拟开采 10 年地下水量均衡分析表Table4Equilibriumanalysisofgroundwateramountinsimulatedminingfor10years均衡项流入流出量/万m3d1流入流出量/万m3年1百分比/%补给项孔隙水降水入渗量6.802482.0083.64河流入渗量0.28102.203.44灌溉入渗量0.85310.2510.46侧向流入量0.2073.002.46小计8.132967.45100岩溶水降雨入渗量0.58211.706.61河道渗漏量2.55930.7529.04侧向径流量0.65237.257.40越流补给量5.001825.0056.95小计8.783204.70100排泄项孔隙水农业开采量1.50547.5023.08越流排泄量5.001825.0076.92小计6.502372.50100岩溶水水源地集中开采量3.501277.5041.18分散开采量5.001825.0058.82小计8.503102.50100孔隙水补排差1.63594.95/岩溶水补排差0.28102.20/024 km预测地下水位模拟区范围城市乡镇河流断裂带图18旧县水源地优化开采至 2030 年地下水流场预测图Fig.18PredictionofgroundwaterflowfieldfromoptimizedexploitationofwaterresourcesinJiuxiancountytill2030952中国岩溶2023年的特征，水位动态变化属降水入渗开采型。水位年际动态呈现出 2 个阶段特征。19902003 年岩溶地下水位基本呈现波动下降；2003 年限采后岩溶地下水位抬升后基本达到稳定。（3）旧县水源地防止发生岩溶塌陷的临界水位为 108m，处于岩溶含水层顶板以上 2m。（4）通过地下水管理模型确定出临界水位，模拟区岩溶水最大可开采量为 8.28.5 万 m3d1，其中，旧县水源地集中开采方式可开采资源量为 3.23.5万 m3d1。参考文献闫佰忠,孙丰博,李晓萌,王玉清,范成博,陈佳琦.气候变化与人类活动对石家庄市藁城区地下水位埋深的影响J.吉林大学学报(地球科学版),2021,51（3）：854-863.YANBaizhong,SUNFengbo,LIXiaomeng,WANGYuqing,FAN Chengbo,CHEN Jiaqi.Impact of climate change andhumanactivitiesongroundwaterdepthofGaochengdistrictinShijiazhuangCityJ.JournalofJilinUniver