某渣场污染物在非饱和岩土介质中的迁移模拟研究

眭素刚 ，徐世光 ，刘文连，范柱国，丁茂斌，韦立德，马宏杰

(1. 昆明理工大学 国土资源工程学院，云南 昆明，650093；2. 中国有色金属工业昆明勘察设计研究院，云南 昆明，650224；3. 云南省地矿局，云南 昆明，653100；4. 中国科学院 武汉岩土力学研究所，湖北 武汉，430071)

对污染物在非饱和土壤中迁移转化规律的研究始于20世纪80年代，学者们在研究非饱和带水分运动的基础上，开始研究污染物在非饱和土壤中的迁移规律。在最初阶段所采用的示踪剂大都是盐类或保守性物质，通过大量的室内及野外土柱物理模拟试验，确定非饱和带垂向一维弥散系数和衰减系数[1-2]。此后，随着研究工作的深入，逐步开始研究重金属在非饱和带的迁移转化规律，考虑土壤液相和固相浓度的分配系数，并借助于Henry，Freundlich和Langmuir的等温吸附模式来表示液相和固相浓度吸附和解吸问题[3-4]。我国对污染物在非饱和土壤中迁移转化的研究已开始重视，如杨金忠等[5]用NaCl作为污染物的示踪剂，采用电导率法测定土柱内的浓度分布；吕家珑等[6]对土壤中磷的运移进行研究；王亚男等[7]对含磷污水淋滤条件下土壤中磷迁移转化进行了模拟试验。我国学者的研究大部分是在确定性模型基础上进行的一维垂向和均质的溶质运移室内研究，这些成果在自然界实际土壤中应用时，有时偏差较大。国内许多研究者应用土壤及地下水溶质运移研究成果，结合尾矿库自身的特点对此进行研究。马腾等[8]在对我国南方某大型尾矿库库区水文地质条件概化的基础上，讨论了不同条件下 U(VI)在地下水中的迁移情况，建立的是二维计算模型，岩土介质采用的是饱和渗流模型。徐晓春等[9]通过对铜陵凤凰山矿林冲尾矿库中重金属元素的空间分布特征及相关土壤、水系沉积物和植物中重金属元素含量变化的研究，探索尾矿重金属元素的空间变化和迁移规律，分析其对周围环境的影响及潜在危害性。张建军[10]以某尾矿坝为研究对象，岩土介质采用饱和渗流模型，采用数值模拟方法建立二维计算模型，并结合现场监测资料，定量分析矿浆中铜离子迁移规律。为此，本文作者应用土壤及地下水溶质运移研究成果，结合尾矿库自身的特点，对尾矿库中重金属元素的空间分布特征及相关土壤、水系沉积物和植物中重金属元素含量变化进行研究，了解尾矿重金属元素的空间变化和迁移规律。目前，通过试验研究和数值模拟等方法，已经取得了一定成果，但其对溶质运移分析都是建立在饱和岩土介质渗流模型基础之上。某渣场使用至今已50 a，库容已接近使用极限，渣场堆渣含有氟和锌等元素，并长期无屏障堆放，其有害物质受大气降水淋滤直接或间接进入渣场下覆黏土层中，对渣场周边地下水可能形成一定的污染。降雨引起该渣场污染物运移是典型的三维溶质迁移问题。通过对研究区内非饱和土壤弥散、渗流特性试验，考虑非饱和渗流过程的溶质迁移数值模拟方法，并运用国际商用软件COMSOL Multiphysics对研究区建立污染物迁移模拟计算模型，运用 COMSOL Multiphysics的标准模型，即无限制蓄水层垂直断面上的稳态渗流以及污染物瞬态传输，模型使用流体流动和溶质传输的多物理场耦合，应用达西方程和对流扩散方程，分析随时间推移污染物在非均质介质中的运移规律[11]，即分析渣场当前污染物传播范围与浓度分布规律，并预测将来污染物传播范围与浓度分布规律，以便为渣场的闭库治理提供科学依据。

1 研究区概况

1.1 水文地质条件

场地内地下水主要类型为坡、残积黏土孔隙水和碳酸盐岩岩溶水。孔隙水在区间处于盆地边缘地带，含水性差，含微量上层滞水，一般多为雨季含水。大气降水是场地内坡、残积黏土孔隙水的主要补给来源。大气降水进入渣场后，在渣体内以垂直运动为主，到达坡、残积土层顶面，由于场地内坡、残积黏土具微透水性，所以，孔隙水在空间上以沿黏土顶面水平运动为主，向下运动为辅，向低基准面径流，或沿沟谷边缘排泄。由于场地内坡、残积黏土层厚为 1.8～29.2 m，一般层厚达4.0 m以上，平均厚度为10.9 m，呈连续分布，直到盆地边缘，即场地内孔隙水垂直下渗受阻，所以，场地内上部坡、残积黏土孔隙水与场地内下部碳酸盐岩岩溶水水力联系较弱，岩溶水补给范围广，埋藏深，场地地下水埋深达120 m以上。通过各种形态岩溶及构造裂隙接受补给并向深部运移、汇集和储存，岩溶水的运动受断裂构造影响，总体上，由北往南运动。

1.2 渣场堆积物性质

渣场主要堆积在沟谷内，沟谷横断面呈“U”字型，沟谷呈西北高南东低。堆积渣呈北西—南东走向，渣场走向长约438 m，最宽处约380 m。钻孔结果表明：渣场堆积物成分主要是熔炼渣①，其次为薄层施工弃土①1。熔炼渣厚度为0.30～36.60 m，整个渣场均有分布，熔炼渣粒径大于 2 mm的质量分数平均为32.3%，小于总质量的50%，但达到总质量的25%，而径粒大于0.075 mm的质量分数平均为85.8%，超过总质量的50%，因此，根据粒径组类的划分标准[12]，场地内堆积熔炼渣普遍为砾砂类(局部地段为粗砂类)；施工弃土一般呈可塑状态，厚度为0.20～5.40 m，在渣场局部分布，但平面上不连续，无规律性。渣场下伏场地主要由黏土②、黏土②1及灰岩③组成。黏土②呈硬塑状态，局部可塑状态，厚度为0.60～29.20 m，整个场地均有分布；黏土②1呈可塑状态，厚度为0.4～2.1 m，在场地局部地段分布，并且都是分布在黏土与灰岩的接触面上。

为了对渣场堆存的堆积物的性质进行鉴别，对渣场熔炼渣①、施工弃土①1、黏土②1及黏土②共4个典型地层取样，并对其浸出毒性试验(共对15个元素及其化合物(Cr+，Hg，F，Ag，CN，Se，Ba，Be，Cd，Cr，Cu，Ni，Pb，Zn和 As)进行鉴别)，判定具有浸出毒性特征的危险物有熔炼渣①和施工弃土①1[13]，其质量浓度超标元素为F-和Zn(表1)。15个元素试验结果中只有F-和Zn质量浓度超过了GB 8978—1996《污水综合排放标准》最高允许排放质量浓度[14]，并且浸出毒性试验质量浓度远远大于现场所测的渣场堆积体淋滤水的质量浓度，所以，该渣场污染物元素及其化合物为F-和Zn。

表1 样品浸出液中危害成分质量浓度Table 1 Concentration limits in sample leaching solution harm ingredients mg/L

由于渣场堆积体以干堆为主，含水量小，蒸发后基本不含水，渣场堆积体对附近土壤、地下水的污染主要通过因大气降水而形成堆积体淋滤水的渗透、吸附、迁移，所以，对渣体淋滤水的试验尤为重要。为取得渣体淋滤水的污染迁移物的真实质量浓度，本研究对渣体淋滤水进行试验。其试验方案如下：在渣体内埋设PVC管，PVC管直接埋深到渣体底部，PVC管底部0.80 m以上全部打上过滤孔，收集渣体在降雨条件形成的淋滤水，PVC管顶部密封，以1 a为试验周期，在降雨后7，15和30 d收集淋滤水进行试验，试验结果(表2)表明：降雨后15 d与30 d淋滤水试验结果相近，说明在降雨15 d后，渣体淋滤水污染物质量浓度达到稳定，比较接近实际渣场堆积体淋滤水的质量浓度。

表2 淋滤水样分析结果Table 2 Analysis result of leaching water sample mg/L

2 非饱和土壤弥散和渗流特性试验

2.1 非饱和土壤弥散特性试验

土壤水动力弥散系数，即土壤溶质扩散弥散系数，是应用数学模型对土壤水盐、农药、重金属等运动进行定量分析的必不可少的参数之一。非饱和土壤水动力弥散系数Dsh(θ)的测定采用水平土柱法，本次试验共选取了黏土②及熔炼渣① 2种主要材料进行试验，试验结果见表3和4。

由表3和4可见：黏土②及熔炼渣① 2种土壤水扩散率D(θ)和非吸附性溶质 Cl-1的水动力弥散系数Dsh(θ)与土壤含水率θ呈显著幂函数关系(图1)，对同一土壤而言，随着含水率的增加，水动力弥散系数Dsh(θ)及土壤水扩散率D(θ)呈幂函数关系递增；在相同含水率条件下，水动力弥散系数Dsh(θ)大于土壤水扩散率D(θ)；在试验含水率范围内，熔炼渣①的D(θ)和Dsh(θ)分别大于黏土②的D(θ)和Dsh(θ)。试验结果反映了2种不同土壤这一溶质运移参数的变化特征，这种特征与土壤水分的运动有密切关系，即在非饱和土壤水分过程中，入渗的速率直接影响溶质的运移状况，在相同含水量条件下，熔炼渣①的溶质运移速率要大于黏土②的溶质运移速率。

表 3 熔炼渣①D(θ)和 Dsh(θ)试验结果Table 3 Test result of slag ① D(θ) and Dsh(θ)

表 4 黏土②D(θ)和 Dsh(θ)试验结果Table 4 Test result of clay ② D(θ) and Dsh(θ)

2.2 非饱和土壤渗流特性试验

土的水分特征曲线定义了压力水头与土壤含水率之间的函数关系，它和水力传导率决定了水分和化学物质在土壤中的运移速度，通过土水特征曲线可获得非饱和土的渗透函数、抗剪强度等有关参数，其可靠性直接影响着模型的预测结果[15-16]。本次非饱和土壤渗流特性试验共选取了黏土②及熔炼渣① 2种材料，试验数据见表5和表6。本次试验选取van Genuchten提出的模型(VG模型) 对试验数据进行拟合(图 2和图3)。

表5 黏土②试样土-水特征曲线测试试验数据Table 5 Test data of clay ② soil-water characteristics curve

图 1 D(θ)、Dsh(θ)与 θ关系曲线Fig.1 Relationship curve of D(θ), Dsh(θ) and θ

表6 熔炼渣①试样土-水特征曲线测试试验数据Table 6 Test data of slag ① soil-water characteristics curve

图2 黏土②试样土-水特征曲线图Fig.2 Clay ② soil-water characteristics curve graph

图3 熔炼渣①试样土-水特征曲线图Fig.3 Slag ① soil-water characteristics curve graph

土的水分特征曲线受多种因素的影响，不同质地的土壤，其水分特征曲线各不相同，差别比较明显。黏粒含量越高，其比表面和可能的结合水质量分数越高，在相同含水率下的吸力越大，在相同吸力条件下储水能力越强，脱水速度越慢[15]。非饱和土中水的渗透主要与吸力相关，当吸力增加时，渗透系数下降。图2与图3所示分别是黏土②与熔炼渣① 2种土的体积含水率与吸力水头的关系比较结果，可以看出：在相同含水率下，黏土②的吸力远大于熔炼渣①，同一吸力水头下黏土②土样的持水性能远远高于熔炼渣①土样的持水性能。

3 渣场污染物迁移模拟计算

为模拟污染物浓度的时空变化规律，并预测污染地下水的瞬时动态与可能影响的范围，对渣场地下环境渗流场三维数值模拟分析和污染物运移场三维数值模拟分析，研究污染物浓度分布规律。本次数值模拟采用包含污染物运移数值模拟算法的国际商用软件Comsol Multiphysics进行[17-18]。渗流场按照饱和-非饱和稳定渗流模拟，污染物运移场按照非稳定溶质运移模拟。

3.1 计算模型

熔炼渣①和施工弃土①1是污染源，施工弃土①1呈不连续分布且含量少，将熔炼渣①和施工弃土①1合并为同一地层建模，设为材料 1；系列地层②以黏土②为主，黏土②1与黏土②渗透系数在同一数量级，将这些地层合并为同一地层建模，设为材料 2；因地下水位埋葬深，近50 m厚度范围内灰岩不受地下水位变动影响，将这个范围内的灰岩作为同一地层建模，设为材料 3；受地下水位影响范围的灰岩作为同一地层建模，设为材料 4；灰岩③以下为渗透系数低的隔水层岩石作为一地层建模，设为材料 5。采用了贯穿渣场的10个剖面建模，模型长为1.840 km，宽为0.89 km，最大高度差为 0.439 km, 生成 14 620个节点67 968个单元(图 4)。

图4 三维计算几何模型(从东南往西北看，单位：m)Fig.4 3-dimension geometric model (seen from southeast to northwest)

3.2 渗流问题初始及边界条件确定

渣场顶部接受降雨垂直入渗，属于第二类的给定通量边界。北边界地下水位以下部分以场地北部静止水位(埋深130 m)对应的空间水位设置水头边界条件，地下水位以上部分按照该处是分水岭可以设置为对称边界(隔水边界)；东边界地下水位以下部分以场地东部静止水位(埋深 120 m)对应的空间水位设置水头边界条件，地下水位以上部分按照该处是分水岭可以设置为对称边界(隔水边界)；南边界地下水位以下部分以场地南部静止水位(埋深 100 m)对应的空间水位设置水头边界条件，地下水位以上部分按照该处是分水岭可以设置为对称边界(隔水边界)；西边界地下水位以下部分以场地西部静止水位(埋深 148 m)对应的空间水位设置水头边界条件，地下水位以上部分按照该处是分水岭可以设置为对称边界(隔水边界)。

对于渣场，考虑到降雨后部分水分很快入渗到渣体下部和黏土交界部位储存起来，上部渣体覆盖层厚而蒸发量少，因此，渣场顶部入渗流量为渣场所在地区年降雨量(1 011.8 mm)，顶部入渗平均速度为3.2×10-8m/s。对于渣场外围地面，考虑到上覆层黏土渗透系数小、入渗率低，存在蒸发量和地表径流等因素，根据经验取顶部入渗流量为场地年降雨量的 5%，顶部入渗平均速度为1.6×10-9m/s。

3.3 污染物弥散问题初始及边界条件确定

渣场顶部设置通量边界，其量值近似设置为(3.2×10-8)·ρ，ρ为表2中CK-3试样污染物的质量浓度。北边界设置为0 mg/L浓度边界，南边界设置为平流通量边界，其余边界设置为0通量边界。材料1中水的初始浓度为x，其余材料初始浓度为0 mg/L。

3.4 计算材料参数

黏土层渗透系数主要参考了室内、外试验结果，并适当考虑了重力耦合作用。计算中材料1和材料2设置为饱和-非饱和渗流特性材料，它们的土水特征曲线如图2和图3所示，相对渗透率与有效饱和度相同，其他材料设置为饱和渗流特性材料(表7和表8)。

表7 渗流计算采用的材料参数表Table 7 Material parameter table used in seepage calculation

表8 弥散试验计算采用的材料参数Table 8 Material parameter table used in dispersion test calculation

3.5 计算结果分析

根据渗流计算结果，不考虑吸附作用，进行污染因子的运移计算，模拟时间分别为30，52，60和100 a。由氟和锌浓度变化分布云图可以看出(图 5)：随着时间推移，具有渣场溶液浓度的溶液慢慢往下、东、西方向迁移，往下迁移到地下水附近，氟和锌弥散速度很大，导致扩散很快，从而再往下位置氟和锌浓度不再因时间推移而增大；100 a内氟质量浓度等于0.001 kg/m3，沿着地表在水平方向向东(西)迁移距离不足170 m(图6和图7)；100 a内Zn质量浓度等于0.000 5 kg/m3，沿着地表在水平方向向东(西)迁移距离不足100 m。从不同剖面对比情况可以看出：虽然具体数量上有不同之处，但是都有共同点；随着时间推移，具有渣场溶液浓度的溶液逐渐往下和东(西)方向迁移，往下迁移到地下水附近，氟和锌弥散速度很大，导致扩散很快，从而再往下位置氟和锌浓度不再因时间推移而增大。100 a内渣场氟迁移到灰岩中岩溶水的浓度小于0.001 kg/m3，即100 a内渣场下灰岩中岩溶水的氟质量浓度满足Ⅱ类水标准[19]；100 a渣场Zn迁移到灰岩中地下水的质量浓度小于0.000 5 kg/m3，即100 a内渣场下灰岩中地下水的Zn质量浓度满足Ⅱ类水标准[19]。

图5 20-20′剖面氟100 a预测质量浓度分布云图Fig.5 Predicted density distribution cloud graph of 20-20′section fluorine at 100 a

表9 渣场周边水样试验结果表Table 9 Test result table of water sample surrounding slag dump

图6 27-27′剖面往西100 m处氟100 a质量浓度分布云图Fig.6 Density distribution cloud graph of 27-27′ section westward 100 m fluorine at 100 a

图7 27-27′剖面往东100 m处氟100 a质量浓度分布云图Fig.7 Density distribution cloud graph of 27-27′ section eastward 100 m fluorine at 100 a

为验证渣场对周边附近地表和地下水现状污染情况，分别在渣场附近取地表和地下水进行水质分析试验(表9)。试验结果表明：渣场周边地表和地下水属于Ⅱ类水范围[19]，即渣库现状下(已经运行52 a)并未对周边地表水和地下水形成污染。这也验证了采用Comsol Multiphysics软件对渣场污染物迁移模拟分析的可靠性。

4 结论

(1) 通过对熔炼渣①及黏土② 2种材料非饱和土壤非吸附性溶质的综合扩散弥散系数Dsh(θ)试验，得出2种土壤水扩散率D(θ)和非吸附性溶质Cl-1的水动力弥散系数Dsh(θ)与土壤含水率θ呈显著幂函数关系，对同一土壤而言，随着含水量的增加，水动力弥散系数Dsh(θ)及土壤水扩散率D(θ)呈幂函数关系递增。试验结果反映了2种不同土壤这一溶质运移参数的变化特征，这种特征与土壤水分的运动有密切关系，即在非饱和土壤水分过程中，入渗的快慢直接影响溶质的运移状况，表现在相同含水量条件下，熔炼渣①溶质运移的速率要大于黏土的溶质运移速率。

(2) 通过对熔炼渣①及黏土② 2种材料非饱和土渗流特性试验，同一吸力水头下黏土②土样的持水能力远远高于熔炼渣①的持水能力。

(3) 采用Comsol Multiphysics软件，对渣场污染物进行迁移模拟分析。不考虑吸附作用，进行污染因子的运移计算，模拟时间分别为30，52，60和100 a。渣场污染物在水平方向污染范围不超过170 m，100 a内渣场污染物氟和锌的迁移污染对周边地下水污染浓度满足二类水标准的。对渣场周边地下水进行取水试验，渣场现状下对周边地下水并未形成污染，这也验证了采用Comsol Multiphysics软件对渣场污染物迁移模拟分析的可靠性。

(4) 基于以上对渣场污染物迁移分析研究，对渣场闭库提出如下对策：① 对渣场坡面及顶面进行防渗层覆盖(可采用两布一膜覆盖)，防止由于雨水渗透熔炼渣而形成淋滤水污染源；在渣场周边设置地下截水墙、排水沟，截断渣场内外浅层孔隙水相互渗流，降低渣场下伏土壤含水量。通过防渗层、截水墙等措施，使进入渣场堆积体及堆积体下土体的水量减低到最小，从而有效降低污染物的迁移速度，防止渣场对周边环境和地下水的污染。② 在防渗层上覆土护坡，防止防渗层老化，并作为绿化种植层种植草皮和灌木，形成较好的水土保持效果。

[1] Yamaguchi T, Moldrup P, Yokosy S. Using breakthrough curves for parameter estimation in the convection dispersion model of solute transport[J]. Soil Science Society of America Journal,1989, 53(6): 1634-1641.

[2] Eshel B, Amos N. Estimating transport parameters in soils by a maximum likelihood approach[J]. Soil Sci Soc Am J, 1987,51(4): 870-875.

[3] Charbeneau R J. Groundwater contaminant transport with adsorption and ion exchange chemistry method of characteristies for the case without dispersion[J]. Water Resour Res, 1987,17(3): 705-713.

[4] Laryean K B, Elrick D E, Robin M J L. Hydrodynmaic dispersion involving adsorption during unsaturated transient water flow in soil[J]. Soil Sci Soc Am, 1982, 46(4): 667-672.

[5] 杨金忠, 叶自桐, 贾维钊, 等. 野外非饱和土壤中溶质运移的试验研究[J]. 水科学进展, 1993, 4(4): 245-252.

YANG Jinzhong, YI Zitong, JIA Weizhao, et al. Field experiment on solute transport in an unsaturated soil[J]. Advances in Water Science, 1993, 4(4): 245-252.

[6] 吕家珑, 张一平, 苏仕平, 等. 土壤磷运移研究[J]. 土壤学报,1999, 36(l): 75-82.

LÜ Jialong, ZHANG Yiping, SHU Shiping, et al. Studies on phosphorus transport in soils[J]. Acta Pedologica Sinica, 1999,36(l): 75-82.

[7] 王亚男, 王红旗, 舒艳. 含磷污水淋滤条件下土壤中磷迁移转化模拟试验[J]. 环境科学学报, 2001, 21(6): 737-741.WANG Yanan, WANG Hongqi, SHU Yan. Simulated study on transferring and transforming of phosphorus in soil leached by waste water[J]. Acta Scientlae Circumstantiae, 2001, 21(6):737-741.

[8] 马腾, 王焰新. U(VI)在浅层地下水系统中迁移的反应-输运耦合模拟: 以我国南方核工业某尾矿库为例[J]. 地球科学:中国地质大学学报, 2000, 25(5): 456-461.

MA Teng, WANG Yanxin. Coupled reaction-transport modeling of migration of uranium(Ⅵ) in shallow groundwater system: A case study of uranium gangue site in southern china[J]. Earth Science: Journal of China University of Geosciece, 2000, 25(5):456-461.

[9] 徐晓春, 王军, 李援, 等. 安徽铜陵林冲尾矿库重金属元素分布与迁移及其环境影响[J]. 岩石矿物学杂志, 2003, 22(4):433-436.

XU Xiaochun, WANG Jun, LI Yuan, et al. The distribution and migration of heavy metal elements of Linchong tailings reservoir in Tongling, Anhui Province, and their environment effects[J].Acta Petrologica Et Mineralogica, 2003, 22(4): 433-436.

[10] 张建军. 某尾矿坝污染物迁移的数值模拟研究[D]. 西安: 西安理工大学水电学院, 2007: 1-54.

ZHANG Jianjun. Numerical simulation of contaminant transport of a certain tailing dam[D]. Xi'an: Institute of Water Resources and Hydro-electric Engineering, Xi'an University of Technology,2007: 1-54.

[11] ZHENG Chunmiao, Wang P P. MT3DMS documentation and user's guide[M]. US: Amry Corps of Engineers, 1998: 200-297.

[12] GB 50021—2001, 岩土工程勘察规范(2009年版)[S].GB 50021—2001, Code for investigation of geotechnical engineering(2009 Edition)[S].

[13] GB 5085.3—2007, 危险废物鉴别标准 浸出毒性鉴别[S].GB 5085.3—2007, Identification standards for hazardous wastes Identification for extraction toxicity[S].

[14] GB 8978—1996, 污水综合排放标准[S].GB 8978—1996, Integrated wastewater discharge standard[S].

[15] 卢应发, 陈高峰, 罗先启, 等. 土-水特征曲线及相关性研究[J]. 岩土力学, 2008, 29(9): 2481-2486.

LU Yingfa, CHEN Gaofeng, LUO Xianqi, et al. Study of soil-water characteristical curve and its influential factors[J].Rock and Soil Mechanics, 2008, 29(9): 2481-2486.

[16] 周葆春, 孔令伟, 陈伟, 等. 荆门膨胀土土-水特征曲线特征参数分析与非饱和抗剪强度预测[J]. 岩石力学与工程学报,2010, 29(5): 1052-1059.

ZHOU Baochun, KONG Lingwei, CHEN Wei, et al. Analysis of characteristic parameters of soil-water characteristic curve(SWCC) and unsaturated shear strength prediction of Jingmen expansive soil[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(5): 1052-1059.

[17] Yates S R. An analytical solution for one-dimensional transport in heterogeneous porous media[J]. Water Resource Res, 1990,26(10): 2331-2338.

[18] Geng Q Z, Irard G G, Ledoux E. Modeling of nitrogen cycle and nitrate transfer in regional hydrogeologic systems[J].Groundwater, 1996, 34(2): 293-304.

[19] GB/T 14848—93, 地下水质量标准[S].GB/T 14848—93, Quality standard for ground water[S].