2.5維氧化還原自然電位數(shù)值模擬

陳曉樂，崔益安，謝 靜，陽 兵，柳建新

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室，長沙 410083；3.中南大學(xué) 教育部有色金屬成礦預(yù)測重點實驗室，長沙 410083)

0 引言

自然電位法是一種傳統(tǒng)的地球物理勘探方法，它利用地表或井中自然電位異常識別地質(zhì)與地球物理目標體，在礦產(chǎn)、油氣勘探及垃圾填埋檢測等領(lǐng)域中都有成功應(yīng)用。自然電位法作為一種被動源法，其野外觀測快捷，成本低廉。同時該方法對多種與污染擴散伴生的滲流運動、氧化還原反應(yīng)以及微生物活動信號較為敏感[1-3]，適用于工程與環(huán)境物探領(lǐng)域，尤其是地下有機污染物檢測及監(jiān)測。Naudet[1-2]和Revil[3]通過自然電位法圈定污染位點的氧化還原前沿(redox fronts);Arora[4]和Gallas[5]使用該方法成功監(jiān)測了垃圾填埋場的泄漏;Mao[6]通過沙箱和現(xiàn)場試驗研究了自然電位層析成像，并將其應(yīng)用于地下水修復(fù)和污染羽流圈定;Abbas[7]分析了通過自然電位數(shù)據(jù)反演得到的一個富有機物污染場地的氧化還原場分布;Cui[8-9]進行了沙箱實驗，嘗試使用自然電位法來監(jiān)測地下污染的擴散和遷移。

有機污染區(qū)域觀測到的地表自然電位異常跟地電池(Geobattery)模型的自然電位很相似[1-2，10]，但直接用地電池模型難以解釋這種自然電位異常，原因在于像垃圾填埋場這類有機污染區(qū)域不存在地電池模型中的電子導(dǎo)體礦物來傳遞電子。因此有研究人員提出可能存在微生物提供電子傳導(dǎo)通道并參與形成地電池模型的情況[1-4，10]。在此猜想下，Revil[3]基于地電池模型提出了生物地電池(Biogeobattery)模型來解釋有機污染區(qū)域的自然電位異常,此后生物地電池模型被廣泛應(yīng)用;Fachin[11]開展了物理模擬實驗，在實驗沙箱中埋設(shè)有機污染物成功地形成了生物地電池;Doherty[12]采用生物地電池模型來解釋在一個廢棄煤氣廠的煤焦油有機污染物引起的自然電位異常;Risgaard-Petersen[13]還將生物地電池模型用于海洋淤泥和海底沉積物中的生物電化學(xué)過程研究。以上研究表明，地電池多源模型能有效表征地下復(fù)雜氧化還原反應(yīng)的電位異常特性?；诘仉姵囟嘣茨Ｐ?，開展地下污染物氧化還原場數(shù)值模擬研究，有助于自然電位觀測數(shù)據(jù)的反演解釋，更好地發(fā)揮自然電位法在環(huán)境領(lǐng)域的應(yīng)用效果。

無限單元用于解決截斷邊界問題，于上世紀70年代提出[14]。無限單元可分為三類：①Astley映射無限單元[15]；②Bettess映射無限單元[16]；③Burnett多極開展無限單元[17]。與有限單元相比，無限單元更合適求解無限域或半無限域問題。肖曉[18]、原源[19]實現(xiàn)了二維有限元-無限元耦合法，并將其成功應(yīng)用于2.5維直流電阻率正反演研究中;湯井田[20]、張林成[21]、Xie[22-23]先后開展了三維有限元-無限元耦合法研究，并將其應(yīng)用于地球物理數(shù)值模擬中。常規(guī)有限元施加截斷邊界條件比較麻煩，而且不利于求解場源靠近邊界的情況。早期二維有限元-無限元耦合法的實現(xiàn)均基于矩形格[18-19]，對于復(fù)雜模型的剖分效果不佳，適應(yīng)性不強。這里基于三角形有限元及應(yīng)用最為廣泛的Astley型單向映射無限元，改進了二維有限元-無限元耦合算法，通過均勻半空間點電源模型驗證該算法的正確性及有效性，并將其應(yīng)用于地下有機污染物降解過程所產(chǎn)生的氧化還原場的數(shù)值模擬研究中，探討自然電位法在環(huán)境污染探查、監(jiān)測領(lǐng)域的可行性。數(shù)值結(jié)果表明，筆者所改進的算法行之有效，而自然電位法在環(huán)境污染領(lǐng)域也有著其獨特的勘探優(yōu)勢。

1 2.5維自然電位邊值問題

首先考慮二維模型的邊值問題，在二維地電模型中，電流源仍為三維源，故需利用空間波數(shù)域的狄拉克函數(shù)特性求解2.5維自然電位問題。2.5維自然電位地電模型滿足的邊值問題為[18-19]：

▽·(σ▽U)-k2σU=-I0δ(A) ∈Ω

(1)

(2)

(3)

其中：σ為介質(zhì)電導(dǎo)率；I0為點源電流強度；δ(A)為與點源位置A有關(guān)的δ函數(shù)；n為邊界外法向單位矢量；k為離散波數(shù)；U為波數(shù)域電位；K0(kr)、K1(kr)分別為零階、一階第二類修正貝塞爾函數(shù)；Ω、Γs、?！薹謩e為經(jīng)傅里葉變換后的研究區(qū)域、地表邊界、無窮遠邊界；r為點源點到外邊界任意點的距離(圖1)。

圖1 參數(shù)示意圖Fig.1 The sketch of parameters

在有限元-無限元耦合法中，無需考慮式(3)，而由無限單元充當(dāng)邊界條件。故有限元-無限元耦合法的邊值問題所對應(yīng)的變分問題為：

(4)

δF(U)=0

(5)

式(4)右端的變分為：

(6)

2 有限元-無限元耦合法

有限元-無限元的空間耦合如圖2所示。在有限元區(qū)域，利用三角形有限元的3個節(jié)點構(gòu)造有限單元形函數(shù)，采用該組形函數(shù)φi近似計算點x=(x,z)的波數(shù)域電位為式(7)。

圖2 有限元-無限元的空間耦合Fig.2 The spacial coupling of finite-infinite elements

(7)

將式(7)代入式(6)有：

(8)

對式(8)中的面積分項做變換有：

(9)

對于無限元區(qū)域，采用的是Astley型單向映射無限單元對其進行剖分[15，18-19]。該無限單元的具體映射過程如圖3所示。

圖3 二維四節(jié)點無限元映射過程Fig.3 The mapping process of 2D infinite elements(a)母單元;(b)子單元

在二維四節(jié)點無限單元中，任意位置的空間坐標可表示為：

(10)

(11)

其中:Mi為坐標映射函數(shù)，其表達式為：

(12)

(13)

(14)

(15)

式中：ξi，ηi為局部坐標。經(jīng)過上述坐標映射之后，電位借助無限單元延伸至無限遠處，并且在無限遠處衰減為零。

無限單元中任意位置的波數(shù)域電位可表示為：

(16)

其中：Ui為各節(jié)點電位；Ni為無限元插值形函數(shù)，其表達式為：

(17)

(18)

(19)

(20)

無限元區(qū)域所滿足的偏微分方程及其變分形式，同樣可分別由式(1)、式(2)、式(4)、式(5)和式(8)表示。故無限元區(qū)域的面積分可表示為式(21)。

(21)

式(9)、式(21)中的Φ分別由有限單元及無限元求取。合并式(9)、式(21)，有：

δUT[(K1+K2)U-F]=0

(22)

其中

(23)

(24)

(25)

式(23)與式(24)的耦合過程可表示為[22-23]：

(26)

3 傅里葉正反變換

這里所采用的波數(shù)及其權(quán)值如表1所示。

由式(22)、式(26)有：

KU=F

(27)

其中:K=K1+K2;U為波數(shù)域電位向量;F為源向量。

在不同波速條件下(表1)，求解方程(27)得到波數(shù)域電位。然后由傅里葉反變換求解得到空間域電位：

(28)

其中,u為空間域電位。經(jīng)以上計算，從而得到二維模型中所有節(jié)點的電位分布。

4 數(shù)值算例

4.1 算法可靠性驗證

采用均勻半空間多源模型進行算法可靠性驗證，模型如圖4所示。模型尺度為40 m×40 m，剖分為40×40個四邊形網(wǎng)格。為具有更好的對比性，在有限元-無限元耦合法中，通過將四邊形網(wǎng)格對半剖分，實現(xiàn)簡單的三角形剖分。點電流源幅值均為1 A，波數(shù)及其權(quán)值選擇如表1所示。數(shù)值模擬結(jié)果(圖5)表明，在相同的網(wǎng)格剖分(節(jié)點分布)條件下，有限元-無限元耦合法的計算精度與具有混合邊界條件的常規(guī)有限單元法相當(dāng)，相對于解析解的均方誤差分別為0.022 7、0.019 0。而實施過程中無需計算邊界積分，計算效率更高，尤其在多源動態(tài)模型中更能進一步體現(xiàn)其優(yōu)勢，驗證了算法的正確性及有效性。

表1 模型計算所采用的波數(shù)k及其權(quán)值g[24]Tab.1 Wave number k and its weight g used in the model calculations

圖4 均勻半空間多源模型Fig.4 Homogeneous half-space multi-source model

圖5 地表電位曲線Fig.5 Surface electric potential curves

4.2 氧化還原地電模型

地下有機污染物的降解可視為動態(tài)過程[25-26]。相關(guān)生物電化學(xué)反應(yīng)所構(gòu)建的地電池模型是隨著時間動態(tài)變化的。對于每個靜態(tài)時刻，其地電模型可由圖7近似表征。隨著氧化還原強度的變化，圖7模型中的電流強度發(fā)生動態(tài)變化。在實際問題中，地下有機污染物的降解進程受到許多因素的影響，如土壤中自由空氣的濃度、有機污染物的種類及其濃度、土壤微生物的類型及其數(shù)量分布等。但結(jié)合微生物的周期性代謝特征，整個降解系統(tǒng)的氧化還原反應(yīng)強度的總體趨勢可簡化為：先緩后急，再緩至峰值，最后再逐漸衰減。根據(jù)該規(guī)律，筆者由圖6所示的假設(shè)曲線近似表征整個降解過程的電流源幅值變化，以定性探究該過程的自然電位異常特征。

地電模型(圖7)可近似模擬微生物降解地下有機污染物時在地下水位附近發(fā)生的氧化還原反應(yīng)所引起的自然電位異常。模型中有限單元區(qū)域尺度為20 m×30 m，采用三角形剖分，除地表外的三個邊界由無限單元作為邊界單元充當(dāng)邊界條件，為四邊形剖分。設(shè)定地下水位位于2 m深度處，上層電阻率為10 Ω·m，下層電阻率為100 Ω·m。假定地下水位附近有兩處位置存在微生物降解有機污染物的現(xiàn)象，其中地下水位上部發(fā)生氧化反應(yīng)，積累負電荷，地下水位下部發(fā)生還原反應(yīng)，積累正電荷，從而形成生物地電池模型。通過離散正負點電流源近似正負電荷的區(qū)域性分布，各點電流源幅值的動態(tài)變化如圖6所示。

圖6 電流強度幅值的動態(tài)變化Fig.6 The dynamic variation of current intensity amplitudes

圖7 生物地電池簡化模型Fig.7 The simplified model of biogeobattery

數(shù)值模擬結(jié)果(圖8、圖9)表明，地下有機污染物降解過程能在地面觀測到自然電位負異常，且通過局部負異常極值點能判斷氧化還原反應(yīng)的分布范圍及場點數(shù)量。地表自然電位曲線也能進一步反映地下有機污染物降解過程的變化情況，與圖6的變化趨勢呈正相關(guān)。

圖8 異常及非異常點自然電位隨電流源的變化Fig.8 The variation of self-potential of abnormal and non-abnormal points with current sources

圖9 地表自然電位異常曲線Fig.9 Surface self-potential curves

5 結(jié)論

改進了原有基于矩形網(wǎng)格的二維有限元-無限元耦合算法，將有限元區(qū)域三角化，更適合模擬復(fù)雜地電模型。首先通過均勻半空間模型驗證了該算法的可靠性，然后將其應(yīng)用于微生物降解地下有機污染物過程所產(chǎn)生的氧化還原場的數(shù)值模擬中，研究地表自然電位對該過程的異常響應(yīng)。得到以下結(jié)論：

1) 基于有限元區(qū)域三角剖分的改進算法能更有效地剖分復(fù)雜地質(zhì)模型，在相同的網(wǎng)格剖分條件下，單次計算速度及計算精度與常規(guī)有限單元法相當(dāng)，且方便處理場源信息及邊界條件，更適用于自然電位等多源動態(tài)模型的數(shù)值模擬。

2) 地質(zhì)微生物參與的地下有機污染物降解過程會在地表呈現(xiàn)自然電位負異常，且可通過地表電位的變化情況評估地下有機污染物的降解進度。

3) 動態(tài)模型的模擬結(jié)果表明地表自然電位異常與地下電流強度呈正相關(guān)。