崔少北， 李曉斌， 楊雙安

(1.河南理工大學(xué) 資源環(huán)境學(xué)院，焦作 454000；2.中原經(jīng)濟(jì)煤層(頁巖)氣河南省協(xié)同創(chuàng)新中心，焦作 454000)

0 引言

2000 年以來，我國工業(yè)化快速發(fā)展，與此同時造成的地下水污染和破壞日益加劇，對生態(tài)環(huán)境也帶來了不利影響。據(jù)水利部2016年月發(fā)布的“地下水動態(tài)月報”顯示，在抽查的18個省市區(qū)地下水抽樣中，地下水水質(zhì)不達(dá)標(biāo)的約占60%[1]。在全國范圍內(nèi)，約有90%的城市地下水遭受到不同程度的污染[2]。城市淺層地下水水質(zhì)受各種污染源綜合作用的影響，包括污水處理廠排污、垃圾填埋場的瀝出液、工廠排污、排水管道滲漏、大氣沉降、建筑垃圾及游樂場排污等[3]。因此，能夠準(zhǔn)確及時地對城市地下水污染情況進(jìn)行有效地監(jiān)測，改善城市生態(tài)系統(tǒng)水環(huán)境質(zhì)量，對生態(tài)建設(shè)、工業(yè)發(fā)展及居民生活有重要意義。

現(xiàn)有地下水監(jiān)測的方法主要采用定期在地下水監(jiān)測井中取樣，用水化學(xué)方法分析水樣中污染物成分和濃度，進(jìn)而確定地下水的污染狀況。顯然，該方法費(fèi)時、費(fèi)力，采樣分析結(jié)果具有滯后性，也不能給出地下污染水體的邊界范圍。因此，尋找一種快速有效的方法，能對城市淺層地下水進(jìn)行長期監(jiān)測已經(jīng)迫在眉睫。地球物理方法主要用于探測地球本體及近地空間的介質(zhì)結(jié)構(gòu)、物質(zhì)組成、形成演化，在工程、環(huán)保、礦產(chǎn)等方面有非常廣泛的應(yīng)用，其中電法勘探在地下水勘測方面具有顯著優(yōu)勢。在利用電阻率差異探測淺層地下水方面，劉宏偉等[4]基于水化學(xué)與地球物理法的萊州灣南岸海(咸)水入侵勘查，結(jié)合水化學(xué)法和地球物理探測方法，探索出EH4法和準(zhǔn)三維高密度電法，較好地指示咸淡水體在空間上的分布，框定海(咸)入侵水界線的大致范圍[4];陳超啟等[5]高采用密度電法對北京市區(qū)周邊某廢棄垃圾填埋場垃圾滲濾液下滲情況進(jìn)行勘查分析，取得較好效果。

城市的各種污染主要造成地下水中的硝酸根離子濃度的增加[6]，在一些工業(yè)發(fā)展水平較高的城市和地區(qū)，地下水中致癌性重金屬嚴(yán)重超標(biāo)[7-8]。硝酸根離子濃度的增加和重金屬滲透會使淺層地下水的電導(dǎo)率顯著提高，地下受污染水體的電學(xué)性能與未污染水體及周圍環(huán)境在電阻率等物理特性方面有顯著差異。地球物理電法的供電電源在地下空間中形成的穩(wěn)定電場，可以圈定出地下水低阻異常區(qū)的范圍，推測地下水體受污染情況，追蹤地下水污染源，為淺層地下水治理提供參考。根據(jù)這一特性，筆者針對城市淺層地下水水質(zhì)預(yù)警，提出通過在潛水層甚至更深的地層中布設(shè)足夠長的供電電極，長期監(jiān)測城市淺層水電阻率變化。由于長供電電極可以深埋地下，提供更多的電流進(jìn)入深部地層，其地電場反映有更多深部信息[9]。因此，針對城市淺層地下水賦存條件，設(shè)計不同地質(zhì)物理模型，分別用點(diǎn)電極與長電極供電電源對模型激勵相同電壓，數(shù)值模擬地下水污染前后的電場分布特征，比較兩種電極在監(jiān)測城市淺層地下水污染狀況方面的差異，論證點(diǎn)電極和長電極電源在監(jiān)測預(yù)警城區(qū)淺層地下水體污染方面的可行性。

1 理論與方法

1.1 理論基礎(chǔ)

常規(guī)直流電法采用不銹鋼或銅電極向地下提供直流電場，利用相同的電極測量地面兩點(diǎn)之間電位差，進(jìn)一步反演獲得地下地質(zhì)體的電學(xué)特性。由于電極長度較短，地面電法勘探中近似為點(diǎn)電流源。然而，采用長電極供電時，其長度不能忽略。假設(shè)大地是一個無窮大的均質(zhì)半空間，地下水體為空間內(nèi)低阻異常區(qū)，分別放置點(diǎn)電極和長電極兩種供電電源，分析其在均質(zhì)空間中電位的分布特征。

在地表A(0,0,0)點(diǎn)處放置電流為I的點(diǎn)電極(如圖1)，則電流密度J、介質(zhì)電導(dǎo)率σ與電位u滿足以下關(guān)系：

J=-σ▽u

(1)

在地下的半空間中做任意閉合曲面Г，曲面包圍的區(qū)域為Ω。根據(jù)通量定律,電流I與電位u滿足微分方程：

圖1 點(diǎn)電極空間電位示意圖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是點(diǎn)電源與邊界?！奚系狞c(diǎn)之間的距離，cos(r,n)是r與外法向n的夾角余弦。在地表處，沿地表法線方向沒有電流，即a=0，則式(3)轉(zhuǎn)化為黎曼(Neumann)邊界條件[13]。令▽·(σ▽u)=-2Iδ(A)，構(gòu)造泛函：

(4)

對式(4)作變分處理，點(diǎn)電極無限半空間區(qū)域的邊值問題與下列變分問題等價:

(5)

對長電極的泛函問題，因為長電極電阻率一般遠(yuǎn)小于周圍地層電阻率，可視長電極為等勢體。將其分解為一系列的電流元dl，每個電流元可以近似為一個獨(dú)立的點(diǎn)電極(如圖2)。在柱坐標(biāo)中，電流元dl在B(r,φ,d)點(diǎn)處的電位可由鏡像法得出。同理，長電極L在半空間中B點(diǎn)的電位與長電極2L在全空間中B點(diǎn)的電位相同。所以沿2L對電流元dl積分，可得長電極在半空間中的電位分布。

由于電位關(guān)于Z軸軸對稱，不考慮角度φ，電流元dl在地下均質(zhì)半空間任意B點(diǎn)處的電位可由下式得出：

(6)

其中：J為長電極L的線電流密度。積分得出表達(dá)式如下：

(7)

圖2所示的半空間中，長電極L是一個線性系統(tǒng)[15]，L的泛函可由式(5)積分給出：

(8)

當(dāng)長電極L穿過不同介質(zhì)時，電極表面電流隨周圍介質(zhì)電導(dǎo)率變化而改變。電流密度分布的解析式可以由長電極表面電位近似導(dǎo)出，不同介質(zhì)的長電極近似為等勢體，由式(1)可得：

(9)

式中:Ji為電極電流密度;σ為長電極所在介質(zhì)電導(dǎo)率;i為長電極穿過介質(zhì)的序數(shù)。當(dāng)長電極穿過的N個不同介質(zhì)時，有以下關(guān)系：

(10)

(11)

在式(10)中:Li為長電極在第i區(qū)域的長度;σi為該區(qū)域的電導(dǎo)率。由以上兩式可得第k段長電極上的電導(dǎo)率：

(12)

1.2 有限元分析

上世紀(jì)50年代，有限元方法在彈性力學(xué)領(lǐng)域大量使用，隨后發(fā)展為地球物理學(xué)中重要的正演數(shù)值模擬方法[16-18]。由于非結(jié)構(gòu)網(wǎng)格沒有規(guī)則的拓?fù)潢P(guān)系，組成網(wǎng)格的單元之間相對獨(dú)立，可以構(gòu)建任何復(fù)雜的幾何模型。

采用多重網(wǎng)格方法將空間域剖分為正四面體空間離散單元，全空間中任一單元K的電導(dǎo)率設(shè)為常數(shù)σ，正四面體網(wǎng)格具有四個節(jié)點(diǎn)，使用線性插值[16]，單元K中電位u由如下公式確定：

(13)

式中：ui為單元K上第i節(jié)點(diǎn)處電位的有限元近似解;φi為對應(yīng)的基函數(shù)或形函數(shù)。將式(13)代入泛函式(8)可得其離散格式：

(14)

式中：

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

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

其中:NE為空間剖分的單元總數(shù);P為全局電源分布項，相對于點(diǎn)電極的電阻率有限元方法;單元矩陣A、B多了沿長電極的線積分。對式(14)取一階變分，并令其等于零，即可得到基于節(jié)點(diǎn)有限元方法的線性方程組：

QU=P

(15)

式中：U代表所有的待求未知量，即所有單元節(jié)點(diǎn)處的電位；Q為全局系數(shù)矩陣，一般為對稱正定稀疏矩陣。各單元內(nèi)部電位均可利用單元網(wǎng)絡(luò)節(jié)點(diǎn)處的電位公式(13)插值得到。

2 模型設(shè)計

根據(jù)常見城市淺層地下水的賦存形態(tài)，設(shè)計圓餅狀、條帶狀、球狀三種模型。圓餅狀模型模擬賦存于在地表有廣泛分布的第四系松散沉積物空隙中的潛水層，其形態(tài)受基巖表面特征與其上部自由水面控制。潛水層在近地表沿地面平行展布，層厚度遠(yuǎn)小于延展半徑。條帶狀是常見的地下水賦存形態(tài)，在層間錯動形成的破碎裂隙帶中或者巖溶形成的地下暗河中，地下水常呈條帶狀分布。在巖溶地區(qū)或一些構(gòu)造復(fù)雜的地區(qū)，地下賦水空間發(fā)育成洞室狀，長寬高數(shù)值相近，以球狀模型模擬此地下水形態(tài)。

圖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 球狀(點(diǎn)狀)地下水模型IIIFig.6 Spherical groundwater model

如圖3所示，圓餅狀模型I的基本形態(tài)。中心藍(lán)色區(qū)域為地下水污染區(qū)，其電導(dǎo)率σ1取值位0.5 S/m，灰色區(qū)域為未污染潛水層，電導(dǎo)率σ2取值0.1 S/m，上層覆土壤和基底隔水層電導(dǎo)率σ3取值0.01 S/m。布置五條測線：SI-1、SI-2、SI-3、SI-4、SI-5。地表水平測線SI-1的起點(diǎn)坐標(biāo)為E(20,0,0)，終點(diǎn)為F(140,0,0)，4條垂直測線距原點(diǎn)O的距離為：40 m、70 m、100 m、130 m,具體參數(shù)參見表1。

條帶狀是常見的地下水賦存空間，圖5模擬地下暗河的地下水分布特性，中心藍(lán)色區(qū)域為地下水污染區(qū)。在模型II中,地下暗河中心坐標(biāo)為(0,0,-30)，點(diǎn)電極與長電極位于坐標(biāo)原點(diǎn)O(0,0,0)，長電極L穿過地下暗河。在模型上布置四條測線：SII-1、SII-2、SII-3、SII-4。測線SII-1平行地下暗河走向，端點(diǎn)坐標(biāo)為(50,0,0)、(550,0,0)。測線SII-2垂直暗河走向，端點(diǎn)坐標(biāo)為(250,200,0)、(250,-200,0)。兩條垂向測線SII-3、SII-4，地表出漏點(diǎn)坐標(biāo)為(250,0,0)(250,50,0),具體參數(shù)參見表2。

模型III模擬賦水空間發(fā)育成洞室狀的分布特性。在A(30,0,0)點(diǎn)放置兩種電極，同樣布置四條測線：SIII-1、SIII-2、SIII-3、SIII-4。水平測線SIII-1、SIII-2互相垂直，坐標(biāo)為SIII-1:(10,0,0)、(-80,0,0),SIII-2:(0,80,0)、(0,-80,0)。垂向測線SIII-3、SIII-4地表出漏點(diǎn)坐標(biāo)為(0,0,0)(-80,0,0),具體參數(shù)參見表3。

3 數(shù)值計算分析

3.1 圓餅狀地下水模型

對圓餅狀地下水模型中點(diǎn)O點(diǎn)處分別放置點(diǎn)電極和長電極，施加200 V的直流電場，無窮遠(yuǎn)處電位為零。經(jīng)過計算，可得地下水污染前后地表測線SI-1上電位變化特征(圖7)，分析圖7可知,點(diǎn)電極附近地表電位數(shù)值較大，隨距離r的增加電位值急劇衰減，在距供電電極較遠(yuǎn)處，電位曲線趨于平緩。

表1 圓餅狀地下水模型參數(shù)值

表2 條帶狀地下水模型參數(shù)值

表3 球狀地下水模型參數(shù)值

圖7 點(diǎn)電極地表測線SI-1電位曲線圖Fig.7 Point electrode surface survey line SI-1 potential curve

圖8 點(diǎn)電極XZ軸面污染中心區(qū)電位等值線圖(有污染)Fig.8 Point electrode XZ axial pollution center area potential contour map(with pollution)

被污染水體對應(yīng)的地表電位數(shù)值較污染前變小、電位下降加快。這是由于地下水體被污染后，電導(dǎo)率增大，吸收更多電流，地表電位隨之降低。

在模型I中XZ軸面上，作地下水體被污染后的等值線圖(圖8)。圖上由外及內(nèi)分別作：0.2 V、0.25 V、0.3 V、1 V，五個電位等值線。以0.2 V等值線所在等值面旋轉(zhuǎn)3600可得模型IV(圖9)。

同理可得XZ軸切面上兩種供電模式在地下水被污染前后的電位分布(圖10、圖11、圖12)。在點(diǎn)電極供電條件下，電位0.23 V等值線分布于污染水體周圍空間，受水體電導(dǎo)率變化影響較小，空間電位的擾動主要出現(xiàn)在污染區(qū)內(nèi)部。地下水被污染后，0.3 V電位等值線分布區(qū)域減小，污染區(qū)內(nèi)電位值處在0.23 V～0.3 V之間，污染水體接更近等勢體。在長電極供電條件下，XZ軸切面上的電位值遠(yuǎn)大于點(diǎn)電極。地下水被污染后，電場擾動范圍擴(kuò)大，電導(dǎo)率增高使得高電位分布區(qū)明顯擴(kuò)大。沿圖中五條測線繪制電位曲線，作地下水被污染前后各測線上電位變化曲線(圖13～圖16)。

圖9 點(diǎn)電極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 點(diǎn)電極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 點(diǎn)電極地表測線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 點(diǎn)電極垂直測線電位曲線圖Fig.15 Point electrode vertical line potential curve

圖16 長電極垂直測線電位曲Fig.16 Long electrode vertical line potential curve

圖17 點(diǎn)電極地表測線電位曲線圖Fig.17 Point electrode surface line potential curve

圖18 長電極地表測線電位曲線圖Fig.18 Long electrode surface line potential curve

圖19 點(diǎn)電極垂直測線電位曲線圖Fig.19 Point electrode vertical line potential curve

圖20 長電極垂直測線電位曲線圖Fig.20 Long electrode vertical line potential curve

在圖13中，點(diǎn)電極供電條件下，污染區(qū)界面以內(nèi)的電位減小，電位曲線下降速率減緩。將80 m～120 m段電位曲線圖放大，兩條電位曲線在污染區(qū)邊界以外基本重合，地下水被污染后地表電位曲線在污染區(qū)邊界處下降速率變大。這是由于在污染區(qū)邊界處，地下水電導(dǎo)率由σ1變?yōu)棣?，電位下降速率變大，電位曲線出現(xiàn)輕微彎折。在垂直方向上，污染區(qū)邊界(測線SI-4)以內(nèi)電位數(shù)值明顯減小，隨深度的增加電位下降速率減緩，潛水層中電位曲線與X軸近似平行。在污染區(qū)邊界以外，SI-4、SI-5兩條測線上的電位分布變化不大，這說明污染水體對遠(yuǎn)處地下空間電位的影響較小，電導(dǎo)率變化引起的空間電位異常主要分布于被污染水體上部地表以及污染水體所在空間層位。

相對而言，長電極供電條件下，水體污染后地表電位數(shù)值變大，電位下降速率減緩，如圖14所示，在污染區(qū)邊界處電位下降速率有明顯改變。在圖16中，潛水層下界面以上空間電位幾乎不隨深度發(fā)生變化，地下水體被污染后的電位遠(yuǎn)大于污染前。污染水體對電位分布的擾動范圍擴(kuò)大，增大的電位區(qū)域延伸至污染區(qū)邊界處(SI-4)以外。

圖21 點(diǎn)電極地表測線電位曲線圖Fig.21 Point electrode surface line potential curve

圖22 長電極地表測線電位曲線圖Fig.22 Long electrode surface line potential curve

圖23 點(diǎn)電極垂直測線電位曲線圖Fig.23 Point electrode vertical line potential curve

圖24 長電極垂直測線電位曲線圖Fig.24 Long electrode vertical line potential curve

綜合以上分析，對于圓餅狀地下含水層，點(diǎn)電極電場受地下水污染狀況影響較小，電位變化不明顯。長電極穿過地下含水層，地下水體被污染后，電導(dǎo)率變大，電位下降速率減緩，對應(yīng)的地表電位值大幅增加。所以，長電極在尋找地下水以及地下水體污染監(jiān)測方面具有顯著優(yōu)勢。

3.2 條帶狀地下水模型

對模型II模擬計算，繪制電位曲線圖(圖17～圖20)。對比圖17～圖20分析，長電極供電條件下電位數(shù)值遠(yuǎn)大于點(diǎn)電極。兩種電極地表測線SII-1上的電位隨r(測點(diǎn)與原點(diǎn)的距離)的增大逐步降低。地下水體被污染后，點(diǎn)電極附近地表電位變小，電位下降速率減緩，隨r的增大電位又逐漸高于污染前電位，在賦水區(qū)邊界處，電位曲線下降速率加快，邊界以外兩條電位曲線趨向于重合。測線SII-2關(guān)于SII-1對稱，地下電位關(guān)于含水空間走向?qū)ΨQ分布。在垂直測線上，穿過含水層的測線SII-3上的電位曲線，在含水層下界面處曲線出現(xiàn)拐點(diǎn)，下界面以上空間中電位變化較小。在右下角放大圖中可以發(fā)現(xiàn)，當(dāng)測線深度超過下界面時，電位數(shù)值快速下降。

3.3 球狀

在模型III中，點(diǎn)電極和長電極的位置偏離球心，長電極插入含水層。對模型III進(jìn)行模擬計算，可以得到球狀儲水空間地下水電位分布特征。可以發(fā)現(xiàn)，長電極電位值大于點(diǎn)電極。點(diǎn)電極供電條件下，測線SIII-1上的兩條電位曲線，差異微小。SIII-2測線上的電位在球心上部地表位置出現(xiàn)極值，污染前的電位大于污染后，這是由于地下水體對電場的擾動在球心處最大。在長電極供電條件下，地表兩條測線上電位曲線接近重合。測線SIII-1與長電極共面，隨著測點(diǎn)與電極距離的增大，電位曲線逐漸降低。SIII-2測線垂直SIII-1，電位曲線呈對稱狀，在原點(diǎn)處取得極值。觀察圖23，穿過模型III的測線SIII-3上，電位曲線有兩個拐點(diǎn)，球型模型所在空間的電位曲線下降平緩，上下兩端測線上電位變化速率較大。這是由于球形儲水空間中水體電導(dǎo)率較大，電位變化速率較小。測線SIII-4在球心等深度線上的電位曲線出現(xiàn)拐點(diǎn)，這是由于球形儲水空間體積較小，其對空間電場的擾動出現(xiàn)在球心以上空間。

4 結(jié)論

筆者結(jié)合常見的地下水形態(tài)，設(shè)計了三種地下水模型。分別用長電極和點(diǎn)電極提供穩(wěn)定電壓，分析地表電位的差異，比較兩種電極在監(jiān)測地下水污染情況方面的優(yōu)劣。通過分析三個模型電位分布特征，得出以下結(jié)論：

1)長電極和點(diǎn)電極通過提供穩(wěn)定的電場，都可以監(jiān)測地下水污染情況。由于長電極可以將更多電流傳入深部地層，其電場帶有更多信息，對地下水電導(dǎo)率變化更加敏感。長電極形成的電場強(qiáng)度更大，電位差異更加明顯，具有較高的信噪比。點(diǎn)電極供電電流受地表覆土阻隔，對地下水電導(dǎo)率變化不敏感。

2)對于不同類型的地下水模型，其電場分布特征不同。地下水體被污染后，水體對電場的擾動受賦存形態(tài)影響。在對地下水污染情況進(jìn)行監(jiān)測時，合理布置供電電極與測線的位置，可以有效獲取地下水污染情況。

3)根據(jù)水文地質(zhì)特征具體分析，考慮到地下水深度和滲流速度的影響，結(jié)合長電極和點(diǎn)電極電源的電法設(shè)計城市淺層地下水污染狀況監(jiān)測預(yù)計系統(tǒng)，可以長期有效的監(jiān)測地下水污染狀況。