水動力和介質(zhì)非均質(zhì)性影響下DNAPLs運移數(shù)值模擬

韓智穎，王錦國，崔孜銘，楊蘊

(河海大學(xué)地球科學(xué)與工程學(xué)院，南京 211100)

隨著經(jīng)濟的飛速發(fā)展，社會對石油化工產(chǎn)品需求量日益增加，導(dǎo)致每年大量的揮發(fā)性液態(tài)有機化合物從地表泄漏并滲透進入地下，對土壤環(huán)境和地下水資源造成嚴重污染。這些液態(tài)有機化合物一般與水互不相溶，而是以非水相液體(Non-aqueous Phase Liquids，NAPLs)的方式與水相、氣相共存于土壤-地下水系統(tǒng)中組成復(fù)雜的多相流運移問題。根據(jù)NAPLs密度大小與水相密度比較，通?？煞譃檩p非水相液體(Light Non-aqueous Phase Liquids，LNAPLs)和重非水相液體(Dense Non-aqueous Phase Liquids，DNAPLs)[1]。

相較于LNAPLs而言，DNAPLs具有更低的粘滯性和界面張力，其運移規(guī)律更為復(fù)雜多變；在均勻介質(zhì)中，DNPALs在非飽和層中的遷移過程與LNAPLs相似，主要受重力和毛細力的共同影響, 同時產(chǎn)生水平向和豎直向的位移；自然條件下重力的影響遠遠大于毛細力，污染物分布形態(tài)總體表現(xiàn)為向下運移；但由于DNAPLs相對于水相密度較大，到達飽和層中后，地下水動力占主導(dǎo)作用，LNAPLs基本沿水平方向運移[2-3]；而DNAPLs則滲過潛水面繼續(xù)運移，并在地下水動力和重力的協(xié)同作用下最終集聚在飽和含水層底部。所以研究DNAPLs在地下介質(zhì)中的運移分布規(guī)律等具有重要的意義和價值。

地下水動力、泄露速率、降雨入滲、多孔介質(zhì)的非均質(zhì)性等是影響DNAPLs遷移和空間分布的主要外部因素，其中地下水動力和介質(zhì)的非均質(zhì)性影響DNAPL運移規(guī)律更為顯著[4-9]。Erning等[4]采用TOUGH系列軟件TMVOC模塊研究了地下水流速變化對DNAPL入滲擴散行為的影響。研究表明，DNAPL體對應(yīng)用的地下水流速的反應(yīng)比場景的幾何設(shè)置更敏感。即使在其位置、尺寸和結(jié)構(gòu)的最低流速下，DNAPL也受到影響。施小清等[7,8]和Yang等[5]研究了介質(zhì)非均質(zhì)性、水相飽和度、DNAPL泄露速率等對DNAPL運移分布的影響，結(jié)果表明介質(zhì)非均質(zhì)性是污染物達到最大飽和度并集聚的主要原因。

非均質(zhì)性是實際天然土層的空間結(jié)構(gòu)最明顯的特征之一，野外的天然地下介質(zhì)通常存在大孔隙、多裂隙發(fā)育、軟弱夾層、低滲透性透鏡體等結(jié)構(gòu)，有時也夾雜著斷裂等地質(zhì)構(gòu)造，從而造成地下多孔介質(zhì)具有高度復(fù)雜的非均質(zhì)性[10]。且由于自然環(huán)境的水力梯度一般橫跨范圍較大且可能有人為因素的影響，例如抽水井附近地區(qū)水力梯度一般大于正常水力梯度，因此本文采用不同水力梯度和不同位置和大小的透鏡體模擬地下水動力和介質(zhì)非均質(zhì)性以及該兩種因素協(xié)同作用對DNAPLs運移規(guī)律和分布特征等影響。

1 多相流運移數(shù)學(xué)模型

TMVOC是水、土壤氣體和揮發(fā)性有機化合物(VOCs)多元混合物在非均質(zhì)多孔介質(zhì)中三相非等溫流動的數(shù)值模擬軟件，采用FORTRAN語言編寫，是勞倫斯伯克利國家實驗室開發(fā)的TOUGH2通用模擬程序的擴展[11]。TMVOC軟件中的多相流計算體系包含水相、NAPL相以及不可壓縮氣體(NCGs)三部分，多相流在遷移的過程各組分之間會發(fā)生擴散并相互轉(zhuǎn)化，并遵循多相、多組分計算體系中流體和熱流的質(zhì)量和能量守衡方程如下：

(1)

式中，q表示質(zhì)量或熱流；q表示匯和源；vn表示單位體積;τn為流體穿過的面積；Mk為組份k在單位體積內(nèi)的質(zhì)量；Fk為組份k進入單位體積內(nèi)的通量；qk為k組份在單位體積內(nèi)的源匯項；n為單位體積向外的法向量矢量，向內(nèi)指向vn。

采用積分有限差分方法進行空間和時間的離散化，所有計算網(wǎng)格塊的相關(guān)主要熱力學(xué)變量用牛頓-拉夫遜迭代法求解，并假設(shè)三相中每相的運動遵循引入相對滲透概念的Darcy定律，如下：

(2)

式中，β表示相位;kγβ為相的相對滲透率(介于0和1之間);pβ=p+pcβ是相中的流體壓力;pcβ是參考相(通常是氣相)中的壓力和毛細管壓力的總和(毛細管壓力為負)。

每個相在壓力和重力的作用下流動，其中包括相之間的相對滲透率和毛細管壓力的影響。本文選用相對滲透率Stone模型和毛管壓力Parker模型，表達式如下：

相對滲透率Stone模型：

(3)

(4)

(5)

式中，kwr、kgr、knr分別為水相、氣相以及NAPL相的相對滲透率;Sw、Sg、Sn分別為水相、氣相以及NAPL相的飽和度(Sw+Sg+Sn=1);Swr、Sgr、Snr分別為水相、氣相以及NAPL相的殘余飽和度;n為擬合參數(shù)。

毛管壓力Parker模型：

(6)

(7)

(8)

(9)

式中，Sm為液相的殘余飽和度;Pcgn為NAPL相-氣相之間的毛管壓力;Pcgw為水相-NAPL相之間的毛管壓力;Swe和Sle分別為水相和液相的有效飽和度;αng、αnw分別為NAPL相-氣相、水相-NAPL相之間的進氣壓力的倒數(shù);n為擬合參數(shù)[12]。

2 理想模型概況與基本參數(shù)

本文采用TMVOC模擬不同水力梯度和不同位置和大小的透鏡體影響下DNAPLs運移規(guī)律和分布特征。算例概化如下：理想模型為二維尺度，長(X)為100 m、寬(Y)為1 m和高(Z)為15 m的長方體(圖1)，共分為680(17×1×40)個網(wǎng)格，其水平方向和垂向的網(wǎng)格剖分如下(表1、表2)。

圖1 概念模型示意圖

表1 X、Y方向網(wǎng)格剖分

表2 Z方向網(wǎng)格剖分

模型的頂部設(shè)定為具有降雨入滲的大氣邊界，為通量邊界；底部以及前后兩個面(XZ剖面)為零通量邊界，左右邊界(YZ剖面)為定水頭邊界，左邊界水頭高于右邊界。模型在水-氣兩相系統(tǒng)中處于平衡狀態(tài)，水相飽和度為1.0，模型溫度恒為20 ℃，不考慮熱物理相關(guān)過程運算。氯苯作為DNAPLs代表，以點源方式在模型頂部以下0.5 m處均勻注入，泄露速率為 1×10-5kg/s，模擬其在連續(xù)注入3 a 內(nèi)的運移和分布結(jié)果。將地質(zhì)體整體概化為粘土巖相，用SAND表示；大氣邊界概化為ATOMS巖相；低滲透性層狀透鏡體模型概化為粘土巖相，用IMP表示。模型的主要參數(shù)值見表3。

表3 模型參數(shù)的設(shè)置

3 結(jié)果與分析

3.1 水力梯度

首先考慮不同水力梯度對DNAPL運移分布的影響，水力梯度分別設(shè)置為0、0.01和0.02，分別對應(yīng)小、中、大3種流速，并選取0.5 a、1 a、2 a 三個時間節(jié)點XVOCM(液相VOC的總質(zhì)量分數(shù))Y=0橫剖面圖作為實驗結(jié)果展示(圖2)。

圖2 不同水力梯度下0.5 a、1 a、2 a XVOCM Y=0橫剖面圖

模擬結(jié)果表明：i=0時，氯苯在重力的影響下垂向遷移，污染范圍基本呈對稱結(jié)構(gòu)；i=0.01時，隨著水力梯度的增加，在0～0.5 a內(nèi)，垂向入滲占據(jù)絕對的主導(dǎo)地位，氯苯沿注入方向快速向下運移。這主要是由于增大的水流流速減小了氯苯驅(qū)替水相的阻力[15]。隨著時間的增加，污染羽受水力梯度的影響逐漸向右運移；在1 a時，污染羽垂向運移的深度和向下游側(cè)運移的長度基本一致；i=0.02時，當(dāng)接近潛水面時，由于水力梯度較大，氯苯便開始隨水流的方向進行彌散。隨著時間的增加，氯苯基本向右運移，污染羽向下運移的深度和向右運移的長度約為1∶5[13]。

3.2 非均質(zhì)性

大部分的學(xué)者主要通過數(shù)值模擬一般采用滲透率隨機場實現(xiàn)介質(zhì)的非均質(zhì)性[5]或室內(nèi)實驗一般用具有不同性質(zhì)或位置大小的夾層或透鏡體來實現(xiàn)介質(zhì)的層狀非均質(zhì)性[16]。本模擬水力梯度設(shè)置為0，在飽水帶設(shè)置不同位置和大小的層狀低滲透性透鏡體，并選取6月、12月、18月、24月、30月、36月時間節(jié)點XVOCM(液相VOC的總質(zhì)量分數(shù))Y=0橫剖面圖作為實驗結(jié)果展示(圖3)。為了定量分析透鏡體上方飽和度的變化，分別選擇網(wǎng)格251、249、254作為3個監(jiān)測單元，監(jiān)測透鏡體上方飽和度在36個月內(nèi)的變化(圖4(a))。

圖3 飽水層透鏡體9～24月XVOCM Y=0橫剖面圖

模擬結(jié)果表明：結(jié)合圖3和圖4(a)看出經(jīng)過6個月左右后，DNAPL運移到透鏡體位置，DNAPL開始在透鏡體正上方集聚，并形成小型污染池；在24個月后，DNAPL飽和度達到一定時，但并未完全滲透穿過透鏡體，便開始向左右兩側(cè)水平運移，并在重力影響下，形成明顯的污染指?？梢钥闯鐾哥R體對于DNAPL在飽水帶的垂向滲透不僅具有明顯的阻滯作用，并加速DNAPL水平運移，使其運移速度明顯大于垂向速度；繞過透鏡體后向下運移，圖4(b)下游側(cè)入滲深度明顯大于上游側(cè)，并在上下游兩側(cè)形成較為明顯的污染指。透鏡體主要發(fā)揮的是阻滯作用，明顯影響了DNAPL運移的分布、結(jié)構(gòu)特征[14]。

圖4 飽水層透鏡體對DNAPL運移速度和分布范圍的影響

3.3 不同水力梯度下透鏡體對DNAPL運移的影響

通過改變透鏡體在模型中的大小和分布位置(圖5)，并選取21月、24月、27月3個時間節(jié)點XVOCM(液相VOC的總質(zhì)量分數(shù))Y=0橫剖面圖作為實驗結(jié)果展示(圖5)。探討水力梯度和介質(zhì)非均質(zhì)性協(xié)同作用對DNAPL運移分布規(guī)律、結(jié)構(gòu)特征的影響。

圖5 不同位置透鏡體21～27月XVOCM Y=0橫剖面圖

模擬結(jié)果表明：

(1) 有透鏡體存在時，在相同的時間內(nèi)污染物的運移明顯遲緩，27月時由于透鏡體的存在，DNAPL污染物比無透鏡體時少垂向運移了近7 m，證明透鏡體對DNAPL污染物的運移具有明顯的阻滯作用。在21月時，DNAPL污染物已經(jīng)在透鏡體上積蓄形成污染池，飽和度逐漸增加；當(dāng)飽和度增加到一定時，DNAPL開始繞過透鏡體垂向運移，形成了明顯的污染指；當(dāng)設(shè)置了從左至右水力梯度后，上游污染指運移速度明顯大于下游側(cè)，且垂向運移的速度也明顯加快。

(2) 當(dāng)水力梯度和透鏡體同時存在時，水力梯度的存在使得介質(zhì)非均質(zhì)性影響進一步增強，進一步增大了優(yōu)勢通道的出現(xiàn)概率；造成DNAPL結(jié)構(gòu)幾何效應(yīng)的原因主要有兩個方面：一是流動水的輸運效應(yīng)，將DNAPL推向下游方向；二是地下水流速越高，傳質(zhì)速率越高，溶解滲濾路徑的速度比DNAPL池快[11]。在無水力梯度的均質(zhì)介質(zhì)中 DNAPL污染物分布的范圍基本以對稱結(jié)構(gòu)分布，當(dāng)?shù)叵滤畡恿Υ嬖跁r增強了地下多孔介質(zhì)的非均質(zhì)性，導(dǎo)致 DNAPL遷移路徑的空間變異性和不規(guī)則性增強產(chǎn)生DNAPL污染物在透鏡體的影響下更容易產(chǎn)生蓄積和復(fù)雜的繞流，形成形態(tài)多變的污染池和污染指交替組合。

4 結(jié)論

本文采用TMVOC軟件通過改變水力梯度和透鏡體位置和大小模擬水力梯度、介質(zhì)不均勻性以及該兩種因素共同對DNAPL污染物運移的具體影響。主要結(jié)論歸納如下：

(1) 水力梯度對DNAPL在入滲過程中和入滲后的行為有顯著影響。隨著水力梯度的增加，其滲流路徑以及稀釋和置換的DNAPL池逐漸朝水流方向傾斜。對于地下水力梯度較高的現(xiàn)場，重力的影響會被削弱。

(2) 在透鏡體的影響下，DNAPL在低滲透性透鏡體上方滯留集聚并加強了其水平遷移的速度逐漸形成污染池，污染的范圍顯著增大；污染池達到一定飽和度時，DNAPL污染物繞過透鏡體，朝著上下游兩側(cè)形成的優(yōu)勢通道繼續(xù)入滲。

(3) 水力梯度和透鏡體同時存在時，水力梯度的存在使得介質(zhì)非均質(zhì)性影響進一步增強，進一步增大了優(yōu)勢通道的出現(xiàn)概率；地下多孔介質(zhì)的非均質(zhì)性增強造成 DNAPL運移路徑的不規(guī)則性增強，使其在飽水含水層產(chǎn)生集聚和繞流，形成了多種形態(tài)的污染指和污染池交替組合。