2.5维氧化还原自然电位数值模拟

2021-02-02陈晓乐崔益安柳建新

物探化探计算技术 2021年1期

陈晓乐，崔益安，谢静，阳兵，柳建新

(1.中南大学地球科学与信息物理学院，长沙 410083；2.中南大学有色资源与地质灾害探查湖南省重点实验室，长沙 410083；3.中南大学教育部有色金属成矿预测重点实验室，长沙 410083)

0 引言

自然电位法是一种传统的地球物理勘探方法，它利用地表或井中自然电位异常识别地质与地球物理目标体，在矿产、油气勘探及垃圾填埋检测等领域中都有成功应用。自然电位法作为一种被动源法，其野外观测快捷，成本低廉。同时该方法对多种与污染扩散伴生的渗流运动、氧化还原反应以及微生物活动信号较为敏感[1-3]，适用于工程与环境物探领域，尤其是地下有机污染物检测及监测。Naudet[1-2]和Revil[3]通过自然电位法圈定污染位点的氧化还原前沿(redox fronts);Arora[4]和Gallas[5]使用该方法成功监测了垃圾填埋场的泄漏;Mao[6]通过沙箱和现场试验研究了自然电位层析成像，并将其应用于地下水修复和污染羽流圈定;Abbas[7]分析了通过自然电位数据反演得到的一个富有机物污染场地的氧化还原场分布;Cui[8-9]进行了沙箱实验，尝试使用自然电位法来监测地下污染的扩散和迁移。

有机污染区域观测到的地表自然电位异常跟地电池(Geobattery)模型的自然电位很相似[1-2，10]，但直接用地电池模型难以解释这种自然电位异常，原因在于像垃圾填埋场这类有机污染区域不存在地电池模型中的电子导体矿物来传递电子。因此有研究人员提出可能存在微生物提供电子传导通道并参与形成地电池模型的情况[1-4，10]。在此猜想下，Revil[3]基于地电池模型提出了生物地电池(Biogeobattery)模型来解释有机污染区域的自然电位异常,此后生物地电池模型被广泛应用;Fachin[11]开展了物理模拟实验，在实验沙箱中埋设有机污染物成功地形成了生物地电池;Doherty[12]采用生物地电池模型来解释在一个废弃煤气厂的煤焦油有机污染物引起的自然电位异常;Risgaard-Petersen[13]还将生物地电池模型用于海洋淤泥和海底沉积物中的生物电化学过程研究。以上研究表明，地电池多源模型能有效表征地下复杂氧化还原反应的电位异常特性。基于地电池多源模型，开展地下污染物氧化还原场数值模拟研究，有助于自然电位观测数据的反演解释，更好地发挥自然电位法在环境领域的应用效果。

无限单元用于解决截断边界问题，于上世纪70年代提出[14]。无限单元可分为三类：①Astley映射无限单元[15]；②Bettess映射无限单元[16]；③Burnett多极开展无限单元[17]。与有限单元相比，无限单元更合适求解无限域或半无限域问题。肖晓[18]、原源[19]实现了二维有限元-无限元耦合法，并将其成功应用于2.5维直流电阻率正反演研究中;汤井田[20]、张林成[21]、Xie[22-23]先后开展了三维有限元-无限元耦合法研究，并将其应用于地球物理数值模拟中。常规有限元施加截断边界条件比较麻烦，而且不利于求解场源靠近边界的情况。早期二维有限元-无限元耦合法的实现均基于矩形格[18-19]，对于复杂模型的剖分效果不佳，适应性不强。这里基于三角形有限元及应用最为广泛的Astley型单向映射无限元，改进了二维有限元-无限元耦合算法，通过均匀半空间点电源模型验证该算法的正确性及有效性，并将其应用于地下有机污染物降解过程所产生的氧化还原场的数值模拟研究中，探讨自然电位法在环境污染探查、监测领域的可行性。数值结果表明，笔者所改进的算法行之有效，而自然电位法在环境污染领域也有着其独特的勘探优势。

1 2.5维自然电位边值问题

首先考虑二维模型的边值问题，在二维地电模型中，电流源仍为三维源，故需利用空间波数域的狄拉克函数特性求解2.5维自然电位问题。2.5维自然电位地电模型满足的边值问题为[18-19]：

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

(1)

(2)

(3)

其中：σ为介质电导率；I0为点源电流强度；δ(A)为与点源位置A有关的δ函数；n为边界外法向单位矢量；k为离散波数；U为波数域电位；K0(kr)、K1(kr)分别为零阶、一阶第二类修正贝塞尔函数；Ω、Γs、Γ∞分别为经傅里叶变换后的研究区域、地表边界、无穷远边界；r为点源点到外边界任意点的距离(图1)。

图1 参数示意图Fig.1 The sketch of parameters

在有限元-无限元耦合法中，无需考虑式(3)，而由无限单元充当边界条件。故有限元-无限元耦合法的边值问题所对应的变分问题为：

(4)

δF(U)=0

(5)

式(4)右端的变分为：

(6)

2 有限元-无限元耦合法

有限元-无限元的空间耦合如图2所示。在有限元区域，利用三角形有限元的3个节点构造有限单元形函数，采用该组形函数φi近似计算点x=(x,z)的波数域电位为式(7)。

图2 有限元-无限元的空间耦合Fig.2 The spacial coupling of finite-infinite elements

(7)

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

(8)

对式(8)中的面积分项做变换有：

(9)

对于无限元区域，采用的是Astley型单向映射无限单元对其进行剖分[15，18-19]。该无限单元的具体映射过程如图3所示。

图3 二维四节点无限元映射过程Fig.3 The mapping process of 2D infinite elements(a)母单元;(b)子单元

在二维四节点无限单元中，任意位置的空间坐标可表示为：

(10)

(11)

其中:Mi为坐标映射函数，其表达式为：

(12)

(13)

(14)

(15)

式中：ξi，ηi为局部坐标。经过上述坐标映射之后，电位借助无限单元延伸至无限远处，并且在无限远处衰减为零。

无限单元中任意位置的波数域电位可表示为：

(16)

其中：Ui为各节点电位；Ni为无限元插值形函数，其表达式为：

(17)

(18)

(19)

(20)

无限元区域所满足的偏微分方程及其变分形式，同样可分别由式(1)、式(2)、式(4)、式(5)和式(8)表示。故无限元区域的面积分可表示为式(21)。

(21)

式(9)、式(21)中的Φ分别由有限单元及无限元求取。合并式(9)、式(21)，有：

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

(22)

其中

(23)

(24)

(25)

式(23)与式(24)的耦合过程可表示为[22-23]：

(26)

3 傅里叶正反变换

这里所采用的波数及其权值如表1所示。

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

KU=F

(27)

其中:K=K1+K2;U为波数域电位向量;F为源向量。

在不同波速条件下(表1)，求解方程(27)得到波数域电位。然后由傅里叶反变换求解得到空间域电位：

(28)

其中,u为空间域电位。经以上计算，从而得到二维模型中所有节点的电位分布。

4 数值算例

4.1 算法可靠性验证

采用均匀半空间多源模型进行算法可靠性验证，模型如图4所示。模型尺度为40 m×40 m，剖分为40×40个四边形网格。为具有更好的对比性，在有限元-无限元耦合法中，通过将四边形网格对半剖分，实现简单的三角形剖分。点电流源幅值均为1 A，波数及其权值选择如表1所示。数值模拟结果(图5)表明，在相同的网格剖分(节点分布)条件下，有限元-无限元耦合法的计算精度与具有混合边界条件的常规有限单元法相当，相对于解析解的均方误差分别为0.022 7、0.019 0。而实施过程中无需计算边界积分，计算效率更高，尤其在多源动态模型中更能进一步体现其优势，验证了算法的正确性及有效性。

表1 模型计算所采用的波数k及其权值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 氧化还原地电模型

地下有机污染物的降解可视为动态过程[25-26]。相关生物电化学反应所构建的地电池模型是随着时间动态变化的。对于每个静态时刻，其地电模型可由图7近似表征。随着氧化还原强度的变化，图7模型中的电流强度发生动态变化。在实际问题中，地下有机污染物的降解进程受到许多因素的影响，如土壤中自由空气的浓度、有机污染物的种类及其浓度、土壤微生物的类型及其数量分布等。但结合微生物的周期性代谢特征，整个降解系统的氧化还原反应强度的总体趋势可简化为：先缓后急，再缓至峰值，最后再逐渐衰减。根据该规律，笔者由图6所示的假设曲线近似表征整个降解过程的电流源幅值变化，以定性探究该过程的自然电位异常特征。

地电模型(图7)可近似模拟微生物降解地下有机污染物时在地下水位附近发生的氧化还原反应所引起的自然电位异常。模型中有限单元区域尺度为20 m×30 m，采用三角形剖分，除地表外的三个边界由无限单元作为边界单元充当边界条件，为四边形剖分。设定地下水位位于2 m深度处，上层电阻率为10 Ω·m，下层电阻率为100 Ω·m。假定地下水位附近有两处位置存在微生物降解有机污染物的现象，其中地下水位上部发生氧化反应，积累负电荷，地下水位下部发生还原反应，积累正电荷，从而形成生物地电池模型。通过离散正负点电流源近似正负电荷的区域性分布，各点电流源幅值的动态变化如图6所示。