使用运输模型和地下水年龄模型进行地下水脆弱性评估的现场示范，北京平原_中文

 Abstract

地下水脆弱性通常被认为是含水层对污染的抵抗力，这实质上就是地下水环境脆弱性。同样，地下水系统的可再生能力可以作为地下水资源脆弱性的指标。本文将地下水脆弱性的概念扩展为地下水环境脆弱性和地下水资源脆弱性，并通过考察溶质从地表到地下水位的运移过程（地下水环境脆弱性）和评估地下水系统的可再生能力（地下水资源脆弱性）来评估扩展的地下水总体脆弱性。以北京平原为例，在HYDRUS中建立了一维建模、二维映射模型来模拟溶质在包气带中的运移过程。选择地下水位的运移时间和溶质浓度作为地下水环境脆弱性的指标。同时，利用MODFLOW/MT3DMS中的三维地下水年龄模型来模拟饱和带的年龄分布，并将模拟的年龄作为地下水资源脆弱性的指标。通过交换垂向入渗量和地下水位高程，将两个模型耦合起来，综合考虑各项指标，最终评价出总体脆弱性。结果表明，总体脆弱性可以较好地反映溶质运移特征（包气带地下水环境脆弱性）和地下水系统可更新性（饱和带地下水资源脆弱性）的综合影响。北京平原最脆弱区多分布在农业区，中等脆弱区多分布在抽水点，最不脆弱区多分布在补给区。
Keywords: Beijing Plain · Groundwater vulnerability · Groundwater age · Solute transport

Introduction

地下水脆弱性评估在地下水管理中的重要性日益提高。地下水脆弱性概念最早在20世纪60年代末提出，用来表示含水层抵抗污染的能力（Margat 1968）。脆弱性较高的含水层被认为更容易受到污染。通常，地下水脆弱性可以使用指标体系来评估，例如DRASTIC（Aller等，1987）、SINTACS（Civita和Maio，1997）和AVI（Stempvoort等，1993）。这些指标体系的每个大写字母表示地下水脆弱性的影响因子。例如，在DRASTIC中，这些因素包括地下水深度（D）、补给速率（R）、含水层类型（A）、土壤类型（S）、地形（T）、包气带影响（I）和电导率（C）。指标体系因其应用方便而广受欢迎。然而，过程建模方法也经常被采用并在案例研究中作为替代选择进行讨论（例如，Neukum 和 Azzam 2009；Yu 等人 2010）。

除了污染之外，地下水还因过度开发而面临问题。特别是在中国这样的发展中国家，快速发展的经济需要大量的水资源。因此，在中国的一些地方，地方政府正面临地下水过度开发的严重局面。因此，地下水系统的可持续性受到研究人员的更多关注（例如，Cao 等人 2013）。更新速度或循环周期通常用于表征地下水的可更新性（例如，Li 等人 2008）。更新速度越高或循环周期越短，表明地下水的可更新性越好。近年来，地下水年龄常被用作地下水可更新性的另一个指标（例如，Gleeson等，2010；Hoque和Burgess，2012；Scanlon等，2012；Wang，2013；Wei等，2014）。地下水年龄定义为地下水补给到含水层后的渡越时间。其基本思想是，含水层越靠近补给源，地下水年龄越年轻，更新能力越好。

由于地下水脆弱性的原始概念是含水层对污染的抵抗力（地下水环境脆弱性），为了考虑地下水可更新性，本研究将地下水脆弱性的概念扩展为包括地下水环境脆弱性和地下水资源脆弱性。以北京平原为案例研究区，进行扩展的总体脆弱性评估。在HYDRUS中建立一维建模、二维绘图模型模拟包气带溶质运移过程，选取运移时间和地下水位溶质浓度作为地下水环境脆弱性指标，同时采用MODFLOW/MT3DMS中的三维地下水年龄模型模拟包气带年龄分布，并将模拟年龄作为地下水资源脆弱性的指标。通过交换垂直入渗量和地下水位高程将两个模型耦合在一起，综合考虑各项指标最终评价出总体脆弱性。结果表明，总体脆弱性可以精细地反映溶质运移特征和地下水系统可更新性的综合影响。北京平原最脆弱区多位于农业区，中等脆弱区多位于抽水点，最不脆弱区多位于补给区。总体地下水脆弱性图可以为当地地下水管理提供有效建议。

Introduction to Beijing Plain

北京是世界上发展最快的地区之一，近几十年来，北京人口增长率几乎保持在50万/年，如今人口已超过2000万，GDP增长率近几年保持在10%以上。庞大的人口和快速发展的经济需要大量的水资源。然而，北京平原位于华北平原的北缘（图1），作为我国干旱半干旱地区的一部分，北京目前面临着严重的水资源短缺问题。由于大部分水资源需求依靠地下水开采，地下水年开采量已达到极限。北京平原北、西被燕山、西山环绕，东、南与广阔的华北平原直接相连，地表和含水层倾向由西北向东南倾斜，过去曾有河流由西北流向东南。然而，在过去的几十年里，由于过度开发，几乎所有的河流都变得昙花一现。之前有很多研究关注北京平原的地下水系统（例如，Sun et al. 2012; Yang et al. 2012; Dong et al. 2013; Zhai et al. 2013）。

北京市位于北京平原的中部。随着中国经济的发展，这座城市自 21 世纪初以来得到了极大的扩张。北京周边地区大多是农业区。农民为城市居民提供食物和蔬菜。因此，农业活动对地下水质量的影响也令人担忧。因此，在本研究中，必须同时考虑地下水环境的脆弱性和地下水资源的可再生性。研究区在北部和西部边界沿华北平原边缘切割，但在东部和南部边界沿北京区的行政边界切割。

Methodologies

为了评估地下水的整体脆弱性，我们建立了一个地下水在包气带和饱和带的耦合模型。HYDRUS-1D（Simunek 等人，2009 年）用于模拟包气带的流动和溶质运移过程，而 MODFLOW 和 MT3DMS 用于模拟饱和带的地下水年龄分布。我们通过交换这两个模型的垂直补给率和地下水位高程，构建了一个接口将这两个模型耦合在一起（图 2）。

Flow and Transport Modeling Through Vadose Zone

地下水流经包气带的概念模型如图3所示。降水从地表下渗，水溶解并将肥料或农药从地下带入包气带。一般来说，流动和输送过程主要由土壤的水力性质决定。然而，在农村地区，植物根系吸收水分，必然会影响地下水流动和溶质输送。地下水从地面到地下水位的渗透可以看作是简单的，可以用一维理查德方程来描述。

因此，将研究区域划分为127×114个方格（宽度为1公里），并为每个方格构建一个HYDRUS-1D模型，以模拟每个方格的垂直流动和输送。因此，该一维建模、二维映射模型的输出给出了不同位置溶质输送特性的差异。

对于包气带的一维流动，一维理查德方程为：

其中 xx 为距地表的深度，hh 为含水量，hh 为水头，KK 为土壤水力传导率，SS 为汇/源。在包气带中，hh 和 KK 可以用 hh 表示。有几种模型描述它们之间的关系，例如 van Genuchten (1980)、Durner (1994) 和 Brooks and Corey (1964) 的模型。本文采用 van Genuchten 模型 (1980)。此外，中国科学院南京土壤研究所 (Shi et al. 2004) 采用 SOILPAR (Acutis and Donatelli 2003) 预测土壤的水力特性，以土壤质地为依据。

植物根系吸水是唯一的汇。因此，S(x)S(x) 为根系吸水速率。S(x)S(x) 取决于植物根系分布。通常，植物根系被认为是沿深度分布的纺锤形。但本文将其简化为梯形分布b(x)b(x)（式（2）），其中LrL_r为植物根系的最大深度。

然后可以根据公式 (3) 来评估 S(x)：

其中，α 为植物系数，ETp 为潜在蒸腾速率。ETp 是在某一标准气象条件下，某一假设参考作物的蒸腾速率（Irmak and Haman，2007）。ETp 通常用修正的 Penman-Monteith 方程来计算。幸运的是，Duan 等（2004）已经完成了北京平原的 ETp 图。另一个参数 α 取决于植物类型。植物类型主要由北京平原的土地利用类型（居住区、林地和农业区；图 4）决定。因此，对于每种土地利用类型，其植物参数（例如 α、Lr）可以估算为其相关植物类型的平均值（表 1）。

表 1 典型土地利用类型的植物类型、估计根长 (Lr) 和作物系数 (α)。

Land use	Residential area	Woodland	Agricultural land
Crops or plants	Ignored	Apple, avocado, deciduous fruit trees	Wheat, corn, been, Alfalfa, grass
Estimated average Lr	Ignored	2 m	30 cm
Estimated average α	Ignored	30	20

该入渗流模型下边界条件设定为地下水位饱和含水量边界，上边界设定为年平均降水量下的常通量边界。通过求解对流-弥散方程来模拟溶质运移。运移模型上边界被视为一个假想的单位浓度脉冲。因此，地下水位浓度趋于形成一条突破曲线C(t)C(t)，其峰值浓度为CmC_m。C(t)C(t)达到CmC_m的渡越时间取为特征渡越时间TtT_t。CmC_m和TtT_t均作为地下水环境脆弱性的指标。

地下水年龄建模

北京平原抽水含水层多为第四纪砂岩，一般由西北向东南增厚，地下水补给主要来自山体降水和侧向基流。第四纪含水层通常经验性地分为浅层组和深层组，但实际上，两层组之间并没有明显的界限，而且，浅层组或深层组通常进一步分为含水层和弱透水层，两者之间也没有明显的界限（图5）。在建模中，浅层组分为4层（2个含水层和2个弱透水层），深层组分为5层（3个含水层和2个弱透水层）。各组含水层电导率相同，各组弱透水层电导率也相同（图6、7）。

建模区域也被划分为 127 × 114 个方格（宽度为 1 公里）（图 8）。边界条件通常设置为北部和西部的指定通量边界，南部和东部的第三类水头相关边界。根据集水区水平衡评估北部和西部山区的侧向基流。

北京平原上分布着广泛分布的抽水井，其中大部分位于城市周围，用于市政供水。此外，农业区还有分布更广的农业抽水井。农业抽水率被认为是区域分布的，并从降水补给中减去。由于政府自 2000 年代初以来一直实施严格的抽水许可控制，因此在过去 10 年中抽水率几乎保持不变。因此，流量模型被认为是稳定状态。地下水年龄模型在 MT3DMS 中实施，并根据 9 个观测到的 14C 年龄和 17 个观测到的 3H 年龄进行校准。

Results and analysis

Solute transit time, Tt

由于注入的溶质浓度是假设的单位浓度，因此 Cm 也是假设的，并以百分比表示。Cm 的模式也代表了所有输入参数的总体影响，但 Cm 似乎更依赖于土壤特性（图 10）。在大多数地区，Cm 的变化范围约为 20% 到 50%，这表明，如果溶质以脉冲方式注入，地下水位处的溶质浓度通常比其注入浓度小得多。

Peak solute concentration, Cm

由于 Tt 和 Cm 是所有影响参数的综合结果，因此另一个关注点是找出主要影响因素。分析了 Tt、Cm 与输入参数（例如饱和电导率 (Ks)、残余水含量 (hr)、饱和水含量 (hs)、降水量 (P)、潜在蒸散量 (ETp) 和地下水位深度 (D) 等）之间的相关系数（表 2）。没有关键的相关性；然而，相关性表明 Cm 更多地依赖于土壤水动力学特性，例如 hr 和 Ks。

Groundwater ages

模拟的地下水年龄（图 11）根据 9 个观测到的 14C 年龄和 17 个观测到的 3H 年龄（图 12）进行了校准。年龄分布呈现出从补给区到排泄区越来越老的趋势。最年轻的年龄大多位于侧边界，最古老的年龄出现在抽水井分布最多的洼地锥体上。最古老的年龄可达 40,000 多年。这表明抽水对地下水循环产生了很大的影响。该模型是人为校准的。建模与观测之间的相关性看起来不太好；然而，在各种案例研究中都发现地下水年龄建模很困难，“因为即使在对数对数尺度上，年龄数据的拟合度在视觉上也比水头数据的拟合度差”（Sanford 2011）。校准困难主要是由于对含水层水力特性的估计不准确。无论如何，校准仍然显示模拟年龄的总体趋势与观测年龄的趋势一致。

总体脆弱性评估

参照包气带溶质运移模型和饱水带地下水年龄模型，选取Tt和Cm作为地下水环境脆弱性指标，以地下水年龄A作为地下水资源脆弱性指标，由此可评价总体脆弱性，如公式4所示。

但需要注意的是，还有其他方法可以构建整体脆弱性。本文以公式（4）为例。在这种情况下，相对脆弱性VV可以表示为公式（5）（图13）：

其中，min 和 max 分别为最小值和最大值。
V 图同时考虑了地下水环境脆弱性和地下水资源脆弱性。
通常，V 可分为三个等级。
高脆弱性区域主要分布在农业区（例如，V > 70），中等脆弱性区域与地下水年龄分布较为一致（例如，30 < V ≤ 70），低脆弱性区域主要位于西北部的补给区（例如，V ≤ 30）。

Conclusions

本文将地下水脆弱性概念扩展为地下水环境脆弱性和地下水资源脆弱性，通过考察溶质从地表到地下水位的运移过程评价地下水环境脆弱性，通过地下水年龄建模评价地下水资源脆弱性，反映地下水系统的可更新性。

研究区选择北京平原，划分为127×114个方格，每个方格宽度为1km。在HYDRUS-1D中构建一维建模、二维映射运移模型模拟包气带溶质运移，在MODFLOW和MT3DMS中构建三维地下水年龄模型。通过交换地下水位高程和垂直入渗速率将两个模型耦合起来。选取溶质从地表到地下水位的运移时间和地下水位峰值浓度作为地下水环境脆弱性的指标，选取地下水年龄作为地下水资源脆弱性的指标。

最后，综合考虑地下水环境脆弱性和地下水资源脆弱性，利用这三个指标对地下水扩展脆弱性进行综合评价。总体脆弱性图表明，总体脆弱性可以精细地反映溶质运移特征和地下水系统可更新性的综合影响。北京平原最脆弱区域多分布在农业区，中等脆弱区域多分布在抽水点，最不脆弱区域多分布在补给区。总体地下水脆弱性图可以为当地地下水管理提供有效建议。

致谢
本研究得到国家自然科学基金（Nos. 51409028）和重庆市高校创新团队建设项目（No. KJTD201305）的资助。

References

Acutis M, Donatelli M (2003) SOILPAR 2.00: software to estimate soil hydrological parameters and functions. Eur J Agron 18(3–4):373–377.
Aller L, Bennet T, Lehr JH, Petty RJ, Hackett G (1987) DRASTIC: a standardized system for evaluating groundwater pollution potential using hydrogeological settings. EPA/600/2–87/035. US Environmental Protection Agency, USA.
Brooks RH, Corey AT (1964) Hydraulic properties of porous media, Hydrology Paper No. 3, Colorado State University, Fort Collins.
Cao GL, Zheng CM, Scanlon BR, Liu J, Li WP (2013) Use of flow modeling to assess sustainability of groundwater resources in the North China Plain. Water Resour Res 49(1):159–175.
Civita M, Maio MD (1997) SINTACS: Un sistema parametrico per la valutazione e la cartografia della vulnerabilita` degli acquiferi all’inquinamento. Metodologia and automatizzazione, Pitagora Editrice, Bologna 1997:191.
Dong YH, Xu HZ, Li GM (2013) Wellhead protection area delineation using multiple methods: a case study in Beijing. Environ Earth Sci 70(1):481–488.
Duan YH, Tao S, Li BG (2004) Spatial and temporal variation of reference crop evapotranspiration in Beijing. Chin Agric Meteorol 25(2):22–25 (in Chinese).
Durner W (1994) Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour Res 32(9):211–223.
Gleeson T, VanderSteen J, Sophocleous MA, Taniguchi M, Alley WM, Allen DM, Zhou YX (2010) Groundwater sustainability strategies. Nat Geosci 3:378–379.
Hoque MA, Burgess WG (2012) 14C dating of deep groundwater in the Bengal Aquifer System, Bangladesh: implications for aquifer anisotropy, recharge sources and sustainability. J Hydrol 444–445(11):209–220.
Irmak S, Haman DZ (2007) Evapotranspiration: Potential or Reference? Electronic Data Information Source: ABE 343. Institute of Food and Agricultural Sciences, University of Florida, US.
Li HL, Boufadel MC, Weaver JW (2008) Tide-induced seawater–groundwater circulation in shallow beach aquifers. J Hydrol 352(1–2):211–224.
Margat J (1968) Groundwater vulnerability to contamination. (doc) 68 SGC 198 HYD, BRGM, Orleans, France.
Neukum C, Azzam R (2009) Quantitative assessment of intrinsic groundwater vulnerability to contamination using numerical simulations. Sci Total Environ 408(2):245–254.
Sanford W (2011) Calibration of models using groundwater age. Hydrogeol J 19:13–16.
Scanlon B, Faunt CC, Longuevergne L, Reedy RC, Alley WM, McGuire VL, McMahon PB (2012) Groundwater depletion and sustainability of irrigation in the US High Plains and Central Valley. PNAS 109(24):9320–9325.
Shi XZ, Yu DS, Warner ED, Pan XZ, Petersen GW, Gong ZG, Weindorf DC (2004) Soil Database of 1:1,000,000 digital soil survey and reference system of the Chinese genetic soil classification system. Soil Surv Horizons 45:129–136.
Simunek J, Sejna M, Saito H, Sakai M, Genuchten MTV (2009) The HYDRUS-1D Software Package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably saturated media. Department of Environmental Sciences, University of California Riverside, California.
Stempvoort VD, Evert L, Wassenaar L (1993) Aquifer vulnerability index: a GIS compatible method for groundwater vulnerability mapping. Can Water Res J 18(1):25–37.
Sun F, Shao HB, Wang WQ, Watanabe N, Bilke L, Yang ZS, Huang ZF, Kolditz O (2012) Groundwater deterioration in Nankou-a suburban area of Beijing: data assessment and remediation scenarios. Environ Earth Sci 67(6):1573–1586.
van Genuchten MT (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 44:892–898.
Wang SX (2013) Hydrochemical and isotopic characteristics of groundwater in the Yanqi Basin of Xinjiang province, northwest China. Environ Earth Sci 71(1):427–440.
Wei GX, Chen FH, Ma JZ, Dong Y, Zhu GF, Edmunds WM (2014) Groundwater recharge and evolution of water quality in China’s Jilantai Basin based on hydrogeochemical and isotopic evidence. Environ Earth Sci. doi:10.1007/s12665-014-3257-y.
Yang Y, Li GM, Dong YH, Li M, Yang JQ, Zhou D, Yang ZS, Zheng FD (2012) Influence of south to north water transfer on groundwater dynamic change in Beijing plain. Environ Earth Sci 65(4):1323–1331.
Yu C, Yao YY, Hayes G, Zhang BX, Zheng CM (2010) Quantitative assessment of groundwater vulnerability using index system and transport simulation, Huangshuihe catchment, China. Sci Total Environ 408(24):6108–6116.
Zhai YZ, Wang JS, Teng YG, Zuo R (2013) Hydrogeochemical and isotopic evidence of groundwater evolution and recharge in aquifers in Beijing Plain, China. Environ Earth Sci 69(7):2167–2177.