羊庄盆地地下水硝酸盐污染数值模拟研究

付晓刚 唐仲华 刘彬涛 蔺林林 卜华 闫佰忠 张永礼

摘要：羊庄盆地地下水资源丰富，近年来盆地内地下水硝酸盐污染呈加重趋势。基于Visual Modflow软件建立了盆地地下水水流与溶质运移数学模型，对羊庄断块地下水硝酸盐的迁移与扩散进行了数值模拟，预测了5a和10 a后硝酸盐在地下水中的分布状况以及浓度变化趋势。结果表明：在预测时段，地下水中硝酸盐浓度呈明显上升趋势，模拟区硝酸盐浓度较高的地区主要分布在辘辘井区域，在地下水动力场的影响下，硝酸盐高浓度区域主要向西南方向扩展。氮肥的过量施用和盆地内养殖场粪便的随意堆放是地下水硝酸盐浓度持续增高的主要原因。

关键词：数值模拟：Visual Modflow：硝酸盐污染：羊庄盆地

中图分类号：P641.8

文献标志码：A

doi： 10.3969/j.issn.1000-1379.2018.06.020

羊庄盆地是滕州市乃至枣庄市地下水资源最为丰富的优质饮用水源区[1]，是滕州市和枣庄市新城区的城市供水水源地，地下水资源的开发为当地经济发展和人民生活水平提高做出了重要贡献，但是受人类活动影响，羊庄盆地地下水中硝酸盐污染呈加剧趋势。含硝酸盐的水源对人体危害较大，可导致高铁蛋白症[2]，因此从供水安全角度出发，对羊庄盆地硝酸盐的动态变化规律及影响因素进行全面分析研究显得十分迫切。

影响地下水硝酸盐动态变化的因素复杂，地下水数值模拟可综合水文地质和T农业、生活污染等方面的数据，直接建立其相互间的联系，使得人们能够定量分析和预测含水层中硝酸盐浓度的动态变化规律[3]。因此选用地下水数值模拟软件Visual Modflow建立地下水流和溶质运移耦合模型[4]，对羊庄盆地地下水中硝酸盐的迁移与扩散进行数值模拟，预测含水层中硝酸盐浓度的变化，并对其主要影响因素进行系统分析，以期为地下水硝酸盐污染控制提供依据。

1 研究区概况

羊庄盆地是大型山间构造盆地，地形总体上东北高、西南低，高程为50～ 600 m[1]。盆地四周较封闭，出口狭窄，以岩溶网络为介质，具有输水导水能力和一定调蓄功能，存在基本统一的水动力场，属于典型的封闭式泉排型岩溶水系统。羊庄盆地由辛召断块、山亭断块和羊庄断块3个水文地质子单元组成[1，5-6]，受地形、地貌和地层岩性制约，盆地3个子单元间的水文地质条件差异较大。辛召断块、山亭断块2个子单元位于盆地中上游，多为基岩裸露的低山、丘陵、山间盆地等，地势高，地形起伏大，地下水资源较为贫乏[1]。位于盆地下游的羊庄断块子单元，断裂构造较多，岩溶发育程度高，富水性强，是整个盆地岩溶地下水循环径流和富集储存的主要场所，也是羊庄盆地水文地质单元的主体部分。

由于羊庄盆地内水质取样点集中在羊庄断块，且盆地内的主要水源地集中在该断块，因此根据盆地地下水监测点位、主要污染源以及主要水源地分布情况，最终确定研究范围为羊庄断块水文地质子单元。该子单元具有明确的边界条件和完善的地下水补、径、排、蓄条件，是一个相对独立、较为完整的地下水子系统。大气降水是该断块地下水的主要补给来源，总体径流方向为由东北向西南，最终通过人工开采排泄。新薛河河谷穿越了整个羊庄断块，河水与地下水存在密切的互补互排关系[5，7]。

2 水文地质概念模型

为了建立地下水流及溶质运移数学模型，对实际水文地质条件进行概化，建立水文地质概念模型。根据实际水文地质条件和研究需要，将含水系统概化为非均质各向同性的二维非稳定流[8]。模拟范围为羊庄断块水文地质子单元，总面积约134 km2，其边界条件：西南及南部边界为地表与地下水分水岭相重合的隔水边界：西部边界中段和北段为化石沟隔水断裂带及地下水分水岭组成的零流量边界，南段为通过新薛河和第四系孔隙含水层向盆地外排泄的给定流量边界[9]：东北和北部为以曹王墓断裂为界的侧向径流补给边界[9]；上边界主要接受大气降水、灌溉人渗补给以及新薛河补给和排泄：根据钻孔资料统计，含水层大多集中在埋深220 m以上，其下岩溶发育微弱，因此将220 m埋深的岩层概化为隔水底板，即下边界。

溶质运移模型主要考虑研究区农业施肥和养殖场污水及粪便淋滤水人渗引起的地下水硝酸盐污染，因此将研究区的养殖场处理为污染点源。另外，由于农业生产施肥会导致地下水面源污染，因此在模型中根据已有农田分布情况，处理为净补给溶质浓度。在溶质运移模型中，研究区北边界处理为补给浓度边界，东、西、南边界均处理为零通量边界。此次评价中蒸发浓度按0计算[10]。

3 数学模型始浓度；г1为给定浓度边界；f（x，y，t）为г1边界上的硝酸盐浓度。

4 数值模拟

4.1 地下水流模型

利用Visual Modflow软件进行地下水流数值模拟计算。选择2013年5月26日至2016年9月26日的水位观测资料对模型进行率定。水流模型经过多次运行、反复调整参数，最终使得观测孔处计算水位和实测水位达到了最佳拟合[11]，见图1。

最终研究区降水人渗补给系数划分为6个分区，岩溶含水层渗透系数划分为15个分区，分区情况见图2、图3，研究区1～6分区的大气降水人渗补给系数分别为0.195、0.139、0.188、0.206、0.192、0.100，岩溶含水层水文地质参数见表1。

4.2 溶质运移模型

4.2.1 弥散系数的确定

水动力弥散系数是进行溶质运移计算的最重要参数。为了准确确定含水层的弥散系数，在羊庄盆地进行了地下水弥散试验。试验将七钼酸铵作为示踪剂，投源孔选择曹王墓断裂附近的庄里30号孔，投入七钼酸铵75 kg。在下游南塘、侯庄、羊庄水源地布置了3眼觀测孔进行钼酸根离子监测。弥散试验投源井和观测井位置见图4。

试验采用逐点求参法求取水动力弥散系数。取各观测孔钼酸根离子峰值时刻t1和各观测孔钼酸根离子检出时刻t2及其对应的浓度C1、C2，利用下式计算纵横向水动力弥散系数。式中：D_L为纵向弥散系数：t1为时刻1，即各观测孔钼酸根离子峰值时刻；t2为时刻2，即各观测孔钼酸根离子检出时刻：x为观测孔到投源孔的距离；C1为时刻tl时钼酸根离子含量：C2为时刻t2时钼酸根离子含量：u为地下水实际流速。

根据各观测孔中示踪剂浓度随时间变化的监测数据，绘制各观测孔示踪剂浓度C与监测时间t的关系曲线，见图5。地下水弥散系数计算参数见表2，地下水弥散度计算结果见表3。

4.2.2 模型硝酸盐污染源强确定

研究区畜禽养殖场较多，大多分布在盆地周边灰岩裸露丘陵地区，主要养殖猪、鸡、鸭。这些养殖场污水、粪便就地排放、堆积，均未采取防渗措施，污水及粪便淋滤水人渗导致地下水出现明显的点状、带状硝酸盐污染。从水样分析来看，以往及现在较大的养殖场分布地段或其附近下游一带，地下水均出现了明显的硝酸盐超标现象，养殖场规模越大、就地排泄物越多、持续时间越长，超标越严重。根据多年地下水污染源调查资料，取渗滤液硝酸盐浓度的多年平均值60.0mg/L。鉴于畜禽养殖场常年污染地下水，因此将养殖场概化为连续补给的给定浓度点状补给源赋值于模型中。

据调查，研究区农田主要种植小麦和玉米，化肥主要施用尿素和复合肥，施用量为1 500～1 800 k～（ hm2.a），土壤中硝酸根离子在30 mg/L左右。模型计算中将其向地下水输送的硝酸盐浓度设定为30.0 mg/L，按耕地分布范围以面状补给的形式赋值于溶质运移模型中。

4.2.3 溶质运移模型的识别

依据研究区已有实测水质数据，模拟期选2013年10月到2015年10月，其中以2014年10月硝酸盐实测浓度场识别模型，以2015年10月硝酸盐浓度场验证模型。经过多次参数调整，最终使模型硝酸盐计算浓度等值线与实测浓度等值线达到最佳拟合状态，见图6、图7。最终确定的弥散度分区见图8，弥散度共分为3个区，1～3区取值分别为34.2、12.3、5.6 m。

从图6和图7可以看出，研究区东南部硝酸盐含量较高，并且呈放射状向四周逐步迁移扩散。计算与实测硝酸盐浓度等值线的整体拟合程度良好，说明溶质运移模型中的参数取值是正确的，所建立的模型能较真实预测研究区地下水中硝酸盐浓度场的变化规律[12]。

5 研究区地下水硝酸盐演化趋势预测

用识别验证后的模型，现状各源汇项不变的情况下，预测5a及10 a后的硝酸盐浓度等值线，见图9～图11。预测结果表明：在预测时段，地下水中硝酸盐浓度呈明显上升趋势，模拟区硝酸盐浓度较高的地区主要分布在辘辘井区域，在地下水动力场影响下，硝酸盐高浓度区域主要向西南方向扩展。

6 结语

基于Visual Modflow软件建立了地下水流与溶质运移数值模型，对羊庄断块地下水硝酸盐的迁移与扩散进行了数值模拟，预测了10 a后硝酸盐在地下水中的分布状况以及浓度变化趋势。

由硝酸盐浓度预测结果可知，若不采取措施，羊庄盆地地下水中硝酸盐浓度将呈增大趋势，这主要是当地农民对氮肥的过量施用和盆地内养殖场粪便随意堆放造成的。因此建议采取必要的措施控制农业施肥的面状污染源，并对以往养殖场发生过污染的场地进行土壤和地下水污染修复，以削减残留污染物对地下水可能造成的二次污染。根据预测结果可知，辘辘井区域地下水中的硝酸盐长期存在浓度较高的现象，主要原因是该处以往建有多个规模较大的养殖场，虽然这些对地下水产生污染的养殖场已于2010年前后被清理拆除，但当时未进行必要的土壤修复，而完全靠含水层自身净化修复需要一个很长的时间过程。为此，建议在地下水受到明显污染的地区，尽快采用人工修复方式降低污染组分的含量。

参考文献：

[1] 苏秋实.山东省羊庄水源地地下水资源预测及开发利用研兖[D].济南：山东师范大学，2009：1-2.

[2] 张新民.低碳农业的外部性与市场失灵[J].天津农业科学，2012，18（2）：61-64.

[3] 张翠云，马琳娜，张胜，等.Visual Modflow在石家庄市地下水硝酸盐污染模拟中的应用[J].地球学报，2007，28（6）：561-566.

[4] 孙志浩，王开章，孔凡亮.基于FEFLOW的浅层地下水水环境演化模拟[J].中国农村水利水电，2009（ 10）：5 -7，11.

[5] 卜华，陈占成，张良鹏.饮用水水源地保护区划分研究：以山东羊庄盆地地下水水源地为例[J].地质调查与研究，2008，31（3）：236-241.

[6] 曲英森，郭留福，冯振中，等.山东省枣庄市羊庄盆地水文地质初步勘探报告[R].兖州：山东省地质局第三水文地质工程地质大队.1981：31-33.

[7] 吴爱民，李传谟，文唐章，等.山东省羊庄盆地岩溶水均衡试验报告[R].兖州：山东省地矿局第三水文地质工程地质大队.1988：50-51.

[8] 张俊华，于福荣，李志萍.FEFLOW在黄龙T业园地下水数值模拟中的应用[J].人民黄河，2011，33（6）：69-70.

[9] 张秋兰，崔亚莉，邵景力，等.基于数值模型的山东省羊庄水源地保护区划分浅析[J].水文地质T程地质，2008，35（增刊1）：423 -427.

[10] 郝蓉.巴彥淖尔市狼山山前冲洪积扇水环境质量评价与溶质运移模拟的研究[D].呼和浩特：内蒙古农业大学，2010：48-50.

[11] 冯雯娟，刘国东，甘蓉，等.基于GMS的岩溶区火电厂地下水污染物运移模拟[J].人民黄河，2015，37（3）：89-91.

[12] 崔欣.镇赉县生态环境脆弱区引水灌溉对地下水及环境的影响研究[D].长春：吉林大学，2008：68 - 69.