电法监测城市浅层地下水污染可行性研究

2020-06-10崔少北李晓斌杨双安

物探化探计算技术 2020年2期

崔少北，李晓斌，杨双安

(1.河南理工大学资源环境学院，焦作 454000；2.中原经济煤层(页岩)气河南省协同创新中心，焦作 454000)

0 引言

2000 年以来，我国工业化快速发展，与此同时造成的地下水污染和破坏日益加剧，对生态环境也带来了不利影响。据水利部2016年月发布的“地下水动态月报”显示，在抽查的18个省市区地下水抽样中，地下水水质不达标的约占60%[1]。在全国范围内，约有90%的城市地下水遭受到不同程度的污染[2]。城市浅层地下水水质受各种污染源综合作用的影响，包括污水处理厂排污、垃圾填埋场的沥出液、工厂排污、排水管道渗漏、大气沉降、建筑垃圾及游乐场排污等[3]。因此，能够准确及时地对城市地下水污染情况进行有效地监测，改善城市生态系统水环境质量，对生态建设、工业发展及居民生活有重要意义。

现有地下水监测的方法主要采用定期在地下水监测井中取样，用水化学方法分析水样中污染物成分和浓度，进而确定地下水的污染状况。显然，该方法费时、费力，采样分析结果具有滞后性，也不能给出地下污染水体的边界范围。因此，寻找一种快速有效的方法，能对城市浅层地下水进行长期监测已经迫在眉睫。地球物理方法主要用于探测地球本体及近地空间的介质结构、物质组成、形成演化，在工程、环保、矿产等方面有非常广泛的应用，其中电法勘探在地下水勘测方面具有显著优势。在利用电阻率差异探测浅层地下水方面，刘宏伟等[4]基于水化学与地球物理法的莱州湾南岸海(咸)水入侵勘查，结合水化学法和地球物理探测方法，探索出EH4法和准三维高密度电法，较好地指示咸淡水体在空间上的分布，框定海(咸)入侵水界线的大致范围[4];陈超启等[5]高采用密度电法对北京市区周边某废弃垃圾填埋场垃圾渗滤液下渗情况进行勘查分析，取得较好效果。

城市的各种污染主要造成地下水中的硝酸根离子浓度的增加[6]，在一些工业发展水平较高的城市和地区，地下水中致癌性重金属严重超标[7-8]。硝酸根离子浓度的增加和重金属渗透会使浅层地下水的电导率显著提高，地下受污染水体的电学性能与未污染水体及周围环境在电阻率等物理特性方面有显著差异。地球物理电法的供电电源在地下空间中形成的稳定电场，可以圈定出地下水低阻异常区的范围，推测地下水体受污染情况，追踪地下水污染源，为浅层地下水治理提供参考。根据这一特性，笔者针对城市浅层地下水水质预警，提出通过在潜水层甚至更深的地层中布设足够长的供电电极，长期监测城市浅层水电阻率变化。由于长供电电极可以深埋地下，提供更多的电流进入深部地层，其地电场反映有更多深部信息[9]。因此，针对城市浅层地下水赋存条件，设计不同地质物理模型，分别用点电极与长电极供电电源对模型激励相同电压，数值模拟地下水污染前后的电场分布特征，比较两种电极在监测城市浅层地下水污染状况方面的差异，论证点电极和长电极电源在监测预警城区浅层地下水体污染方面的可行性。

1 理论与方法

1.1 理论基础

常规直流电法采用不锈钢或铜电极向地下提供直流电场，利用相同的电极测量地面两点之间电位差，进一步反演获得地下地质体的电学特性。由于电极长度较短，地面电法勘探中近似为点电流源。然而，采用长电极供电时，其长度不能忽略。假设大地是一个无穷大的均质半空间，地下水体为空间内低阻异常区，分别放置点电极和长电极两种供电电源，分析其在均质空间中电位的分布特征。

在地表A(0,0,0)点处放置电流为I的点电极(如图1)，则电流密度J、介质电导率σ与电位u满足以下关系：

J=-σ▽u

(1)

在地下的半空间中做任意闭合曲面Г，曲面包围的区域为Ω。根据通量定律,电流I与电位u满足微分方程：

图1 点电极空间电位示意图Fig.1 Point electrode space potential diagram

图2 长电极空间电位示意图Fig.2 The long electrode space potential diagram

(2)

可以采用混合边界条件[10-14]，即：

(3)

式中:a=cos(r,n)/r，r是点电源与边界Γ∞上的点之间的距离，cos(r,n)是r与外法向n的夹角余弦。在地表处，沿地表法线方向没有电流，即a=0，则式(3)转化为黎曼(Neumann)边界条件[13]。令▽·(σ▽u)=-2Iδ(A)，构造泛函：

(4)

对式(4)作变分处理，点电极无限半空间区域的边值问题与下列变分问题等价:

(5)

对长电极的泛函问题，因为长电极电阻率一般远小于周围地层电阻率，可视长电极为等势体。将其分解为一系列的电流元dl，每个电流元可以近似为一个独立的点电极(如图2)。在柱坐标中，电流元dl在B(r,φ,d)点处的电位可由镜像法得出。同理，长电极L在半空间中B点的电位与长电极2L在全空间中B点的电位相同。所以沿2L对电流元dl积分，可得长电极在半空间中的电位分布。

由于电位关于Z轴轴对称，不考虑角度φ，电流元dl在地下均质半空间任意B点处的电位可由下式得出：

(6)

其中：J为长电极L的线电流密度。积分得出表达式如下：

(7)

图2所示的半空间中，长电极L是一个线性系统[15]，L的泛函可由式(5)积分给出：

(8)

当长电极L穿过不同介质时，电极表面电流随周围介质电导率变化而改变。电流密度分布的解析式可以由长电极表面电位近似导出，不同介质的长电极近似为等势体，由式(1)可得：

(9)

式中:Ji为电极电流密度;σ为长电极所在介质电导率;i为长电极穿过介质的序数。当长电极穿过的N个不同介质时，有以下关系：

(10)

(11)

在式(10)中:Li为长电极在第i区域的长度;σi为该区域的电导率。由以上两式可得第k段长电极上的电导率：

(12)

1.2 有限元分析

上世纪50年代，有限元方法在弹性力学领域大量使用，随后发展为地球物理学中重要的正演数值模拟方法[16-18]。由于非结构网格没有规则的拓扑关系，组成网格的单元之间相对独立，可以构建任何复杂的几何模型。

采用多重网格方法将空间域剖分为正四面体空间离散单元，全空间中任一单元K的电导率设为常数σ，正四面体网格具有四个节点，使用线性插值[16]，单元K中电位u由如下公式确定：

(13)

式中：ui为单元K上第i节点处电位的有限元近似解;φi为对应的基函数或形函数。将式(13)代入泛函式(8)可得其离散格式：

(14)

式中：

u= (u1,u2,u3,u4)T

φ= (φ1,φ2,φ3,φ4)T

其中:NE为空间剖分的单元总数;P为全局电源分布项，相对于点电极的电阻率有限元方法;单元矩阵A、B多了沿长电极的线积分。对式(14)取一阶变分，并令其等于零，即可得到基于节点有限元方法的线性方程组：

QU=P

(15)

式中：U代表所有的待求未知量，即所有单元节点处的电位；Q为全局系数矩阵，一般为对称正定稀疏矩阵。各单元内部电位均可利用单元网络节点处的电位公式(13)插值得到。

2 模型设计

根据常见城市浅层地下水的赋存形态，设计圆饼状、条带状、球状三种模型。圆饼状模型模拟赋存于在地表有广泛分布的第四系松散沉积物空隙中的潜水层，其形态受基岩表面特征与其上部自由水面控制。潜水层在近地表沿地面平行展布，层厚度远小于延展半径。条带状是常见的地下水赋存形态，在层间错动形成的破碎裂隙带中或者岩溶形成的地下暗河中，地下水常呈条带状分布。在岩溶地区或一些构造复杂的地区，地下赋水空间发育成洞室状，长宽高数值相近，以球状模型模拟此地下水形态。

图3 圆饼状地下水模型IFig.3 Round cake groundwater model I

图4 圆饼状地下水XZ轴切面图Fig.4 XZ axis cutaway view of round cake groundwater

图5 条带状(线状)地下水模型IIFig.5 Strip Groundwater Model II

图6 球状(点状)地下水模型IIIFig.6 Spherical groundwater model

如图3所示，圆饼状模型I的基本形态。中心蓝色区域为地下水污染区，其电导率σ1取值位0.5 S/m，灰色区域为未污染潜水层，电导率σ2取值0.1 S/m，上层覆土壤和基底隔水层电导率σ3取值0.01 S/m。布置五条测线：SI-1、SI-2、SI-3、SI-4、SI-5。地表水平测线SI-1的起点坐标为E(20,0,0)，终点为F(140,0,0)，4条垂直测线距原点O的距离为：40 m、70 m、100 m、130 m,具体参数参见表1。

条带状是常见的地下水赋存空间，图5模拟地下暗河的地下水分布特性，中心蓝色区域为地下水污染区。在模型II中,地下暗河中心坐标为(0,0,-30)，点电极与长电极位于坐标原点O(0,0,0)，长电极L穿过地下暗河。在模型上布置四条测线：SII-1、SII-2、SII-3、SII-4。测线SII-1平行地下暗河走向，端点坐标为(50,0,0)、(550,0,0)。测线SII-2垂直暗河走向，端点坐标为(250,200,0)、(250,-200,0)。两条垂向测线SII-3、SII-4，地表出漏点坐标为(250,0,0)(250,50,0),具体参数参见表2。

模型III模拟赋水空间发育成洞室状的分布特性。在A(30,0,0)点放置两种电极，同样布置四条测线：SIII-1、SIII-2、SIII-3、SIII-4。水平测线SIII-1、SIII-2互相垂直，坐标为SIII-1:(10,0,0)、(-80,0,0),SIII-2:(0,80,0)、(0,-80,0)。垂向测线SIII-3、SIII-4地表出漏点坐标为(0,0,0)(-80,0,0),具体参数参见表3。

3 数值计算分析

3.1 圆饼状地下水模型

对圆饼状地下水模型中点O点处分别放置点电极和长电极，施加200 V的直流电场，无穷远处电位为零。经过计算，可得地下水污染前后地表测线SI-1上电位变化特征(图7)，分析图7可知,点电极附近地表电位数值较大，随距离r的增加电位值急剧衰减，在距供电电极较远处，电位曲线趋于平缓。

表1 圆饼状地下水模型参数值

表2 条带状地下水模型参数值

表3 球状地下水模型参数值

图7 点电极地表测线SI-1电位曲线图Fig.7 Point electrode surface survey line SI-1 potential curve

图8 点电极XZ轴面污染中心区电位等值线图(有污染)Fig.8 Point electrode XZ axial pollution center area potential contour map(with pollution)

被污染水体对应的地表电位数值较污染前变小、电位下降加快。这是由于地下水体被污染后，电导率增大，吸收更多电流，地表电位随之降低。

在模型I中XZ轴面上，作地下水体被污染后的等值线图(图8)。图上由外及内分别作：0.2 V、0.25 V、0.3 V、1 V，五个电位等值线。以0.2 V等值线所在等值面旋转3600可得模型IV(图9)。

同理可得XZ轴切面上两种供电模式在地下水被污染前后的电位分布(图10、图11、图12)。在点电极供电条件下，电位0.23 V等值线分布于污染水体周围空间，受水体电导率变化影响较小，空间电位的扰动主要出现在污染区内部。地下水被污染后，0.3 V电位等值线分布区域减小，污染区内电位值处在0.23 V～0.3 V之间，污染水体接更近等势体。在长电极供电条件下，XZ轴切面上的电位值远大于点电极。地下水被污染后，电场扰动范围扩大，电导率增高使得高电位分布区明显扩大。沿图中五条测线绘制电位曲线，作地下水被污染前后各测线上电位变化曲线(图13～图16)。

图9 点电极XZ轴切面电位分布图(有污染)Fig.9 Point electrode XZ axis cut surface potential distribution map (Pollution)

图10 长电极XZ轴切面电位分布图(有污染)Fig.10 Long electrode XZ axis cut surface potential distribution map (Pollution)

图11 点电极XZ轴切面电位分布图(无污染)Fig.11 Point electrode XZ axis cut surface potential distribution map (non-pollution)

图12 长电极XZ轴切面电位分布图(无污染)Fig.12 Long electrode XZ axis cut surface potential distribution map (non-pollution)

图13 点电极地表测线SI-1电位曲线图Fig.13 Point electrode surface measurement line SI-1 potential curve

图14 长电极地表测线SI-1电位曲线图Fig.14 Long electrode surface measurement line SI-1 potential curve

图15 点电极垂直测线电位曲线图Fig.15 Point electrode vertical line potential curve

图16 长电极垂直测线电位曲Fig.16 Long electrode vertical line potential curve

图17 点电极地表测线电位曲线图Fig.17 Point electrode surface line potential curve

图18 长电极地表测线电位曲线图Fig.18 Long electrode surface line potential curve

图19 点电极垂直测线电位曲线图Fig.19 Point electrode vertical line potential curve

图20 长电极垂直测线电位曲线图Fig.20 Long electrode vertical line potential curve

在图13中，点电极供电条件下，污染区界面以内的电位减小，电位曲线下降速率减缓。将80 m～120 m段电位曲线图放大，两条电位曲线在污染区边界以外基本重合，地下水被污染后地表电位曲线在污染区边界处下降速率变大。这是由于在污染区边界处，地下水电导率由σ1变为σ2，电位下降速率变大，电位曲线出现轻微弯折。在垂直方向上，污染区边界(测线SI-4)以内电位数值明显减小，随深度的增加电位下降速率减缓，潜水层中电位曲线与X轴近似平行。在污染区边界以外，SI-4、SI-5两条测线上的电位分布变化不大，这说明污染水体对远处地下空间电位的影响较小，电导率变化引起的空间电位异常主要分布于被污染水体上部地表以及污染水体所在空间层位。

相对而言，长电极供电条件下，水体污染后地表电位数值变大，电位下降速率减缓，如图14所示，在污染区边界处电位下降速率有明显改变。在图16中，潜水层下界面以上空间电位几乎不随深度发生变化，地下水体被污染后的电位远大于污染前。污染水体对电位分布的扰动范围扩大，增大的电位区域延伸至污染区边界处(SI-4)以外。

图21 点电极地表测线电位曲线图Fig.21 Point electrode surface line potential curve

图22 长电极地表测线电位曲线图Fig.22 Long electrode surface line potential curve

图23 点电极垂直测线电位曲线图Fig.23 Point electrode vertical line potential curve

图24 长电极垂直测线电位曲线图Fig.24 Long electrode vertical line potential curve

综合以上分析，对于圆饼状地下含水层，点电极电场受地下水污染状况影响较小，电位变化不明显。长电极穿过地下含水层，地下水体被污染后，电导率变大，电位下降速率减缓，对应的地表电位值大幅增加。所以，长电极在寻找地下水以及地下水体污染监测方面具有显著优势。

3.2 条带状地下水模型

对模型II模拟计算，绘制电位曲线图(图17～图20)。对比图17～图20分析，长电极供电条件下电位数值远大于点电极。两种电极地表测线SII-1上的电位随r(测点与原点的距离)的增大逐步降低。地下水体被污染后，点电极附近地表电位变小，电位下降速率减缓，随r的增大电位又逐渐高于污染前电位，在赋水区边界处，电位曲线下降速率加快，边界以外两条电位曲线趋向于重合。测线SII-2关于SII-1对称，地下电位关于含水空间走向对称分布。在垂直测线上，穿过含水层的测线SII-3上的电位曲线，在含水层下界面处曲线出现拐点，下界面以上空间中电位变化较小。在右下角放大图中可以发现，当测线深度超过下界面时，电位数值快速下降。

3.3 球状

在模型III中，点电极和长电极的位置偏离球心，长电极插入含水层。对模型III进行模拟计算，可以得到球状储水空间地下水电位分布特征。可以发现，长电极电位值大于点电极。点电极供电条件下，测线SIII-1上的两条电位曲线，差异微小。SIII-2测线上的电位在球心上部地表位置出现极值，污染前的电位大于污染后，这是由于地下水体对电场的扰动在球心处最大。在长电极供电条件下，地表两条测线上电位曲线接近重合。测线SIII-1与长电极共面，随着测点与电极距离的增大，电位曲线逐渐降低。SIII-2测线垂直SIII-1，电位曲线呈对称状，在原点处取得极值。观察图23，穿过模型III的测线SIII-3上，电位曲线有两个拐点，球型模型所在空间的电位曲线下降平缓，上下两端测线上电位变化速率较大。这是由于球形储水空间中水体电导率较大，电位变化速率较小。测线SIII-4在球心等深度线上的电位曲线出现拐点，这是由于球形储水空间体积较小，其对空间电场的扰动出现在球心以上空间。