随着人类社会的快速发展，地表工、农业等活动对浅部地下水质量的影响愈来愈强烈。在中国工农业经济发达的平原区，浅部地下水受到不同程度的污染已经具有普遍性[1-4]。对于分布广泛的大型平原区的孔隙承压水，一般来说由于其上部存在弱透水层，水质通常较好，常常成为区域的主要供水层位，但是由于对其大量超采，使得浅层地下水与深部承压水之间的水力联系不断增强，浅层地下水中的污染物在越流的作用下进入深部承压水，进而对其水质造成极大的威胁[5]。越流含水层系统中由于弱透水层的较大比表面积和低渗透性，地下水在弱透水层中的水-岩相互作用较普通多孔介质更加强烈且持续时间更长，因此对大型平原区承压水进行弱透水层对污染物阻滞作用评价是十分重要的。目前，在平原区区域地下水的敏感性分析和越流含水层系统溶质运移研究方面，已经受到国内外学者们的关注[6-11]。
不同的污染组分在地下水系统中迁移时受到吸附、弥散、水岩作用等多个要素影响程度是不同的，因此，结合特征组分的迁移行为进行越流含水层系统弱透水层溶质运移评价研究也越来越受到关注[12]。越流含水层系统弱透水层溶质运移常考虑的水-岩相互作用有吸附-解吸作用、阳离子交换作用，对于不同的溶质和多孔介质，所需考虑的反应也有所不同[13-15]。由于农业面源氮污染的普遍性，近年来学者们开展了不同水文地质条件下地下水对氮的防污性能评价和氮的溶质运移转化研究[16-22]。虽然前人对防污性评价进行了不同角度的改进，在实践中也取得了较好的效果，但是前人对越流含水层系统弱透水层（一般为粘性土层）的处理方面，仍然存在着不足。承压含水层上部的弱透水层对污染物的阻滞能力对区域承压水防污性能起到决定作用。在当前的承压水防污性评价中，普遍的做法是将弱透水层厚度作为评价因子，赋予相应分值和权重来考虑弱透水层的作用，但弱透水层由于沉积特征、岩性特征和粒度差异，其渗透性差异非常显著，对于同一个区域，不同位置弱透水层渗透能力和对污染物的阻滞能力的差异是显著的，仅靠统一赋予评分和权重的方法，不能客观表征弱透水层在承压水防污性评价中的重要作用。因此，有必要开展弱透水层对污染物的阻滞能力定量评价。
本次研究以通州区为例，针对研究区浅部潜水农业氮污染的普遍性，以及承压水在地区供水中的重要作用，对研究区弱透水层对氮污染阻滞能力进行定量模拟分析。研究结果为该地区的地下水开发和保护提供科学依据。
1 研究区浅部含水层系统概况
首都副中心通州区位于华北平原北部，面积为907 km2（图1），为大陆性季风气候，常年平均气温11.3℃，降水的85%集中在6——9月汛期，多年平均降水量620 mm左右。研究区地形平缓，自西北向东南倾斜，地表为第四系覆盖，第四系地层总厚度在150～500 m之间，在水文地质单元划分上位于潮白河冲洪积扇下游和永定河冲洪积扇的扇缘部位，由于沉积特征的差异，地层岩性和地层结构也具有一定的差异性。
图1 通州区地理位置示意图
Fig.1 The location of Tongzhou area
研究区第四系地下水包括潜水和承压水，二者都具有砂层与粘性土相间的互层结构。潜水含水层组底板埋深为45～60 m，水位标高13 m左右，以中、细、粉砂为主，中间夹弱透水层，层多而薄的砂层间为分布不连续的粘性土层。潜水补给源主要为大气降水、河流入渗、地下侧向径流、灌溉回渗等。浅层承压水，埋藏深度在80～120 m，水位标高8 m左右，共有约4个含水层，累积厚度在35～40 m。区域上，由于潜水水位略高于下部承压水水头，因此水动力条件具有向下越流补给承压水的趋势（图2）。潜水与承压水之间分布有较为连续的粘性土层，弱透水层厚度为2～20 m不等，岩性主要为粉质粘土、粘土、粉土构成[17-23]。
图2 研究区第四系浅部地层剖面
Fig.2 The shallow Quaternary stratigraphic profile
由于工农业活动，研究区潜水普遍受到不同程度的污染，深部承压含水层成为当地生活用水的主要供水水源。由于潜水水位高于承压水位，水质较差的潜水向下越流可能造成下部承压水污染。前期调查发现，NH4+、NO3-为潜水污染的重要组分，因此，定量评价弱透水层对特征组分NH4+、NO3-的阻滞能力具有重要意义[24-25]。
2 典型地层介质中NH4+、NO3-迁移室内试验
在资料搜集整理分析和野外实际调查的基础上[26-28]，为明确野外实际地层中含有的污染组分——NH4+、NO3-本底值，及各类岩性的弥散系数、吸附系数和最大吸附量，进行了相关室内实验：（1）基于构建的浅层地下水系统典型含水层组代表类型，开展NH4+、NO3-淋洗实验，测定NH4+、NO3-在地下水含水层介质中的本底值，其中NO3-浓度介于7.11～19.10 mg/kg，NO2-浓度介于0～0.15 mg/kg，NH4+浓度介于0.03～0.14 mg/kg。（2）采用具有典型代表性的钻孔采集的原状土，通过实验获取了通州区典型地层的孔隙度、渗透系数、弥散度等参数，为数值模拟研究NH4+、NO3-在典型地下水系统模式中的迁移提供基础参数（表1）。（3）为研究NH4+、NO3-在地层中的吸附作用，开展了静态吸附和动态吸附实验。静态吸附实验针对土壤对NH4+、NO3-的吸附能力，确定典型土壤对其静态吸附量和相关参数（表2）由高到低依次为粉土、细砂和中砂；动态吸附实验是研究NH4+、NO3-在土壤垂直方向上的迁移行为（表3），从分配系数上看，吸附能力由高到低依次为粉土、细砂和中砂。通过上述室内实验测得了相关参数，为数值模拟奠定了基础。
表1 实验获取孔隙度、渗透系数、弥散度参数表
Table 1 The parameters of porosity,permeability coefficient and dispersity
岩性中砂中细砂粉土孔隙度n/%36.86 32.77 25.33渗透系数m/d 17.32 4.18 0.2弥散系数m2/d 0.484 1 0.112 7 0.002 436弥散度m 0.024 21 0.011 27 0.004 06
表2 静态吸附实验结果
Table 2 Results of static adsorption experiments
土壤岩性中砂细砂Langmiur吸附方程1 S=0.037 1 C+0.003 1 S=0.082 1 C+0.0014吸附系数K L/kg 0.08 0.02最大吸附量Sm mg/kg 333.33 714.29
表3 动态吸附实验结果
Table 3 Results of dynamic adsorption experiments
岩性中砂细砂粉土土样质量g 746.06 864.76 692.37土柱总吸附量mg 230.44 417.41 418.28最大吸附量Sm mg/kg 308.9 543.02 604.12分配系数Kd L/kg 7.72 13.58 15.1
3 承压水上部弱透水层对NH4+、NO3-的阻滞能力模拟探索研究
3.1 模拟方案
针对研究区弱透水层对特征污染组分NH4+、NO3-的阻滞性能，以区域水文地质单元及富水性为基础，结合钻孔资料获取的地层结构特征、岩性特征等，进行潜水-弱透水层-承压水地下水含水层组中特征污染组分NH4+、NO3-迁移模拟方案构建，共构建15种模拟方案。各个方案的对应的水文地质单元分区见图3，表4。
表4 模拟方案
Table 4 The simulation schemes
方案编号水文地质单元分区岩性水力梯度氨氮浓度mg/L硝酸盐浓度mg/L 1 2 3 I I II细砂粉质粘土中砂中粗砂粉土中粗砂中粗砂粉质粘土中细砂厚度/m 15.80 9.40 13.10 20.60 16.50 17.60 15.70 18.84 14.60 0.24 0.98 1.11 0.11 0.98 1.11 0.06 1.21 7.84 4 5 6 7 8 9 10 11 12 13 14 15 II III III III IⅤIⅤIⅤⅤⅤⅤIⅤIⅤII粉细砂粘土细砂中细砂粉土中细砂中细砂粉土中砂粉细砂粉土细砂中细砂粉土中砂中细砂亚粘夹粉砂中砂中砂夹粉土粉土细砂细砂夹粉土粉质粘土中细砂细砂粉粉土夹细砂细砂中细砂粉质粘土夹粉土中细砂细砂粉土中细砂夹粉土粉细砂夹粉土粉土中细砂粉质粘土粉砂中砂粉土细砂粉土粉土细砂粉质粘土粉土中细砂粉土粉细砂粉土11.30 14.80 9.10 16.60 18.40 24.00 9.70 15.00 12.10 12.70 13.10 11.10 10.50 11.2 18.30 9.40 14.10 4.70 18.30 25.70 8.90 11.30 21.30 12.40 3.80 12.80 10.50 4.00 23.70 3.50 15.80 13.40 10.90 2.70 14.60 19.40 17.30 32.20 8.60 26.40 6.10 14.10 15.20 0.08 0.06 0.16 0.12 0.09 0.07 0.09 0.17 0.16 0.19 0.14 0.04 1.21 1.21 0.42 0.42 0.42 0.79 0.79 0.79 0.79 0.87 0.87 0.87 7.84 7.84 18.76 18.76 18.76 0.005 0.005 0.005 0.005 49.54 49.54 49.54
3.2 弱透水层阻滞性能定量模拟
（1）水文地质概念模型
模型目标层位为浅层地下含水层组，包括潜水含水层、弱透水层和承压含水层。由于模拟的主要目的为分析弱透水层对污染组分NH4+和NO3-的阻滞能力，故采用垂向一维溶质迁移模型。水文概念模型示意图见图3。含水层组设为饱和状态，并将上、下水流的边界条件设为一类边界。含水层组上部的浓度边界设为浓度已知的边界，浓度赋值为各模拟方案的代表分区的地下水NH4+和NO3-的平均浓度值，其下部浓度边界设为零浓度的梯度边界。典型污染物自上向下迁移过程中考虑“三氮”的相互转化作用。由于NO2-作为“三氮”转化的中间产物，具有不稳定性，故弱透水层阻滞性能分析主要针对NH4+和NO3-进行。为了分析污染物到达含水层组不同位置的浓度随时间变化过程，在模型中潜水含水层的底部和承压含水层的顶部位置设置了两个监测点M1，M2（图4）。通过分析污染物到达监测点M1和M2的时间差，可以获得位于承压含水层上部的弱透水层对特征污染组分阻滞能力的定量化信息。
图3 通州区第四系水文地质单元划分[26]
Fig.3 The division of quaternary hydrogeological units in Tongzhou area
图4 研究区水文地质概念模型
Fig.4 The sketch map of the hydrogeological conceptual model
（2）数学模型
整个反应过程主要涉及水流运动模型和溶质运移模型。
水流运动模型：
其中：θ为体积含水量[L3L-3]；t是时间[T]；h是压力水头[L]；c为常数；K是Ⅴan-Geneuchten模型中的水力传导系数[LT-1]。
溶质运移模型：
铵根
硝酸根
式中：下标1，2分别表示铵根和硝酸根；w和s分别表示氨的溶液态和固态；c是氮素浓度[ML-3]；ρb是土壤容重[ML-3]；Dw是土壤中溶质的弥散系数[L2T-1]；q是达西流量[LT-1]；
和
分别是相似一阶反应速率[-]；μw是零阶反应速率[-]；Kd是液相和固相的吸附系数[L3M-1][29]。
（2）模拟软件及模型参数
模拟软件选用目前国际流行的HYDRUS-1D，该软件广泛应用于饱和-非饱和土壤中水分、盐分、污染物运移研究[29]。模型中各个参数的确定过程：研究区内各个介质的水分特征参数根据其颗粒分析数据进行赋值，并运用HYDRUS软件的神经网络预测法算得，包括饱和含水率（θs）、残余含水率（θr）、饱和水力传导系数（Ks）、进气值的倒数（α）、孔径分布系数（n）和孔径连通系数（l）等。介质密度（ρb）根据测试获得，分配系数（Kd）、水动力弥散参数（Dw）由室内实验测得，且在经验值允许范围内进行调节；其他硝化（k1）、亚硝化（k2）和反硝化（k3）系数则根据前人研究的经验值范围内调节。具体参数值见表5。
表5 模型设置参数
Table 5 Model parameters
参数岩性粉质粘土粉土粉细砂细砂中砂中粗砂中细砂Kd cm3/mg 0.11 0.10 0.09 0.04 0.004 0.003 0.008 ρb g/cm3 1.75 1.73 1.70 1.66 1.60 1.58 1.64 Dw （cm2/d）4.32 4.32 4.32 4.32 4.32 4.32 4.32硝化k1（1/a）0.001 0.002 0.003 0.003 0.004 0.004 0.003亚硝化k2（1/a）36.50 36.50 36.50 36.50 36.50 36.50 36.50反硝化k3（1/a）0.70 0.65 0.25 0.28 0.32 0.35 0.30
（3）模拟结果分析
以方案1为例，简要分析硝酸氮和铵氮模拟的结果。方案1的模拟结果见图5，图6，代表了细砂（潜水）-粉质粘土（弱透水层）-中砂（承压水）的模式（表1）中的NH4+、NO3-的迁移过程。由图5显示，NH4+由潜水含水层顶部运移到其底部需14.6年，再穿过弱透水层进入下部承压含水层的顶部需56.2年。由图6显示，NO3-由潜水含水层顶部运移到其底部需0.11年，再穿过弱透水层进入下部承压含水层的顶部需0.22年。根据方案1的弱透水层厚度（9.40 m），可计算出在方案1模型代表区域平均水力梯度条件下，弱透水层单位厚度中铵氮和硝酸氮的迁移速率分别为0.23 m/a和86.36 m/a，弱透水层阻滞能力与之成反比。
图5 NH4+变化曲线
Fig.5 NH4+change curves
图6 NO3-变化曲线
Fig.6 NO3-change curves
平原区地层系统结构复杂，本次针对通州区构建了15种地下水含水层组模式，研究了NH4+、NO3-的时空迁移转化特征。对于北部和东部概化为单层地层结构的区域，模型中忽略了砂层内薄层不连续的弱透水层，因此，模拟结果在这些地区迁移速度可能偏快。对于通州南部和西部含水层为互层结构的区域，潜水含水层主要为细砂、粉砂和粉细砂，弱透水层多为粘土、粉质粘土，下部承压含水层为粉细砂、细砂。实际地层的互层结构在模拟中进行了一定程度的合并处理，实际的地层结构其防污性能要优于模拟结果。15种方案模拟结果统计见表6。
表6 不同方案弱透水层对NH4+、NO3-的阻滞能力统计表
Table 6 The statistical table of the aquitard retard capacity of different schemes on NH4+and NO3-
水文地质单元分区I I II II III III III IV IV IV V V VI VI VII方案编号1 2 3 4 5 6 7 8 9 10 11 12 13 14 15弱透水层对典型污染物单位厚度的阻滞能力（m/a）硝酸根NO3-86.36 58.93 44.19 36.25 115.63 75.00 68.42 73.33 33.64 47.92 27.66 36.49 60.87 109.38 42.42铵根NH4+0.23 0.36 0.08 0.06 0.13 0.24 0.21 0.24 0.11 0.25 0.1 0.17 0.17 0.08 0.09
4 结论
根据表6数据结合水文地质分区图3可见，位于研究区东南部水文地质II区、ⅤI区、ⅤII区的潞城-西集-漷县镇一带由于区域弱透水层介质粘粒含量多、厚度大，其对污染物的阻滞能力强防污性能较好；位于西北部水文地质I区、III区的马驹桥镇-台湖镇西部一带弱透水层介质中粗颗粒含量相对较多，厚度相对较小，弱透水层污染物阻滞能力最差。总体来看，研究区弱透水层污染物阻滞能力西北部最差，自西北向东南呈现出先增强后减弱的态势。
鉴于弱透水层在承压水防污性中的关键作用，本次利用数值模拟的方法定量化了弱透水层对氮污染的阻滞作用，下一步会将其阻滞作用作为承压水防污性评价指标体系的因子，并结合研究区环境和水文地质特征，构建研究区孔隙承压水的防污性评价指标体系。然后应用该评价体系进行承压水防污性评价。该研究的思路和方法，可以为浅部孔隙承压水的防污性能评价提供参考和借鉴。
[1]席北斗，李娟，汪洋，等.京津冀地区地下水污染防治现状问题及科技发展对策[J].环境科学研究，2019，32（1）:1-8.
[2]苗晋杰，肖国强，谢海澜，等.天津滨海新区供水结构及对策分析[J].地质调查与研究，2008，31（1）:12-15.
[3]薛禹群，张幼宽.地下水污染防治在我国水体污染控制与治理中的双重意义[J].环境科学学报，2009，29（3）:474-481.
[4]张新钰，辛宝东，王晓红，等.我国地下水污染研究进展[J].地球与环境，2011，39（3）:415-421.
[5]覃蓉高，曹广祝，仵彦卿.非均质含水层中渗流与溶质运移研究进展[J].地球科学进展，2014，29（1）：30-41.
[6]郭高轩，琚宜文，翟航，等.北京市通州区地下水分层质量评价及水化学特征[J].环境科学，2014，35（6）:2114-2119.
[7]苗晋杰，靳继红，杜东，等.首都副中心及重点区域地下水环境质量评价与问题成因[J].地质调查与研究，2020，43（03）：224-229.
[8]费宇红，张兆吉，郭春艳，等.区域地下水质量评价及影响因素识别方法研究-以华北平原为例[J].地球学报，2014，35（2）：131-138.
[9]张兆吉，费宇红，郭春艳，等.华北平原区域地下水污染评价[J].吉林大学学报（地球科学版），2012，42（5）:1456-1461.
[10]郭高轩，琚宜文，翟航，等.北京市通州区地下水分层质量评价及水化学特征[J].环境科学，2014，（06）:2114-2119.
[11]刘丽雅.浑河傍河区地下水氮污染来源贡献识别[D].中国地质大学（北京），2013.
[12]FIJANI E,NADIRI A A,ASGHARI MOGHADDAM A,et al.Optimization of DRASTIC method by supervised committee machine artificial intelligence to assess groundwater vulnerability for Maragheh-Bonab plain aquifer,Iran[J].Journal of Hydrology，2018，503,89-100.
[13]成建梅，胡进武.饱和水流溶质运移问题数值解法综述[J].水文地质工程地质，2003，（2）:99-106.
[14]宋爽.多环芳烃（菲）在弱透水层中越流迁移及衰减规律研究[D].长春:吉林大学，2011.
[15]闫金忠，梁贵荣，邢述彦.地下水溶质运移中的化学反应[J].太原理工大学学报，1998，29（1）：437-440.
[16]RUPERT M G.Calibration of the DRASTIC ground water vulnerability mapping method[J].Groundwater,2001,39,625-630.
[17]NOLAN B T.Relating nitrogen sources and aquifer susceptibility to nitrate in shallow ground waters of the United States[J].Ground Water，2010，39,290-299.
[18]TESORIERO A J,VOSS F D.Predicting the Probability of Elevated Nitrate Concentrations in the Puget Sound Basin:Implications for Aquifer Susceptibility and Vulnerability[J].Ground Water,1997,(35):1029-1039.
[19]PACHECO FERNANDO A L,FERNANDES S LUIS F.The multivariate statistical structure of DRASTIC model[J].Journal of Hydrology，2013，(476):442-459.
[20]Mclay C D A,Dragten R,Sparling G,et al.Predicting groundwater nitrate concentrations in a region of mixed agricultural land use:a comparison of three approaches[J].Environmental Pollution,2001,(115):191-204.
[21]白耀楠，马震，张竞，等.廊坊北三县工程建设适宜性评价[J].地质调查与研究，2019，42（02）:117-122.
[22]张玉川.北京市通州区工程建设层地质建模与环境质量评价[D].中国地质大学（北京），2009.
[23]黄骁，王进卫，陈刚，等.通州规划新城岩土地基工程能力研究[J].城市地质，2015，（01）:25-30.
[24]蔡向民，栾英波，郭高轩，等.北京平原第四系的三维结构[J].中国地质，2009，（05）:1021-1029.
[25]郭高轩，琚宜文，翟航，等.北京市通州区地下水分层质量评价及水化学特征[J].环境科学，2014，（06）:2114-2119.
[26]北京市水文地质工程地质公司.北京市通县卫星城城区地下水资源勘查报告[D].北京市水文地质工程地质公司，1987.
[27]北京市地质矿产勘查开发局.北京市通州区北京通州规划新城前期区域工程地质勘查报告[D].北京市地质矿产勘查开发局，2014.
[28]北京市地质工程勘察院.北京市通州区农村安全饮水工程第四系地下水资源调查评价报告[D].北京市地质工程勘察院，2007.
[29]KIRKHAM AMES M,SMITH CHRISTOPHER J,DOYLE RICHARD B,et al.Inverse modelling for predicting both water and nitrate movement in a structured-clay soil(Red Ferrosol)[J].Peerj(2019),2167-8359.