飽和砂土中DNAPL污染物遷移過程及數值模擬

肖 鵬,劉漢樂(桂林理工大學,廣西環境污染控制理論與技術重點實驗室,巖溶地區水污染控制與用水安全保障協同創新中心,廣西 桂林 541004)

重非水相液體(DNAPL)主要與電鍍、化工和農藥等工農業有關,常見的DNAPL 污染物有三氯乙烯(TCE)、氯代溶劑等,如TCE 可以通過呼吸和皮膚接觸等方式進入人體,可能導致人體部分器官癌變,對人類的潛在危害較大,由于其低水溶性、弱遷移性與難降解性,TCE 能夠穿透含水層而滯留在含水層底部,形成長期的DNAPL 污染源區[1],因此,對于飽和介質中DNAPL 的遷移規律和分布特征的研究得到了國內外學者的重視[2-3].

近年來,國內外學者開展了大量關于DNAPL 在多孔介質中運移規律的研究,由于野外復雜的地質條件與較高的監測成本,多數研究以室內實驗尺度為主.其中,基于透射光學圖像法進行了一系列多孔介質二維砂槽中DNAPL 的遷移和捕集實驗研究,確立了該類方法在其他槽實驗中應用的潛力[4-5].圖像法作為一種快速、準確的非侵入式監測方法,可有效研究DNAPL 在多孔介質中的遷移機理[6-7].已有研究利用圖像法對DNAPL 在含透鏡體的多孔介質中的遷移規律進行了探討,研究表明DNAPL 運移至透鏡體處時無法克服介質間毛細壓力進入透鏡體中,停止垂向入滲行為并在透鏡體表面累積,后在毛細力影響下開始橫向擴散,繞過透鏡體后繼續向下遷移[8-9].而對于野外環境中DNAPL 遷移規律的研究,現有問題在于野外地質環境條件復雜,目前僅限于個別點位的鉆孔取樣分析,對于三維空間DNAPL污染范圍的評估和預測缺乏有效工具,雖然圖像分析法能有效獲取污染物在二維多孔介質區域中的遷移規律,但在定量分析污染物飽和度分布特征上無法應用于野外的三維環境.因此,為進一步探究DNAPL 在飽和地下水中的污染入滲過程,需結合數值模擬方法來確定入滲過程中DNAPL 的位置與空間分布特征.利用室內實驗的優勢,在可控條件下模擬DNAPL 污染入滲的過程,為多相流數值模擬的驗證提供可靠數據[10-12].

DNAPL 在進入地下水飽和帶后,地下介質孔隙空間處于油-水兩相平衡,此時DNAPL 在含水層中的運移受多因素的影響,如介質的非均質性[13-16]、地下水流速[17-19]、DNAPL 泄漏條件[20]和DNAPL 自身的物理化學性質等[21].因此,在地下環境中DNAPL 入滲行為的研究上,相較于有一定局限性的室內試驗,方便、經濟適用且能分析建立不同地質條件的數值方法被廣泛使用.許多通過數值模擬的方法研究地下環境中DNAPL 的遷移行為,以此來驗證補充室內實驗.施小清等通過多相流數值模擬方法探究在隨機滲透率場和重力影響下的DNAPL 遷移過程,發現DNAPL 污染物的空間質心位置與污染擴散范圍在垂向上的影響程度要高于水平方向[22].Yang 等[23]通過多相流數值模擬研究揭示了地下異質性對DNAPL 的吸附和質量分布具有強烈影響,地下環境中水飽和度的變化與非均質性的增強均會導致DNAPL 入滲深度的增大;在此基礎上,Wu等[24]模擬研究地下水系統中DNAPL 遷移過程的污染羽流演化,其演化過程同樣受到非均質性的強烈影響.綜上,利用數值模擬方探究DNAPL 在飽水基質中的運移與分布規律越來越受到關注.

PetraSim 數值模擬程序是適用于TOUGH 系列代碼的交互式預處理器和后處理器,可模擬多孔和裂縫介質中多組分、多相流體的非等溫流動,包括水、空氣和揮發性有機化學品的三相流.相比TOUGH 系列其它的模塊,PetraSim 的優點是可快速建立模型并運行檢查模型的模擬結果.缺點是對野外復雜情況進行模型構建時,要求調整的參數數據量龐大,且模擬結果多偏向理想情況.對于PetraSim的實際應用情況,國內外已有許多研究.在國內,郭亞敏采用PetraSim 模擬程序建立了在研究區內的DNAPL 三相流運移模型,有效模擬了土壤中的DNAPL 殘留相在重力與水流作用下對巖溶含水層的入滲行為[25].在國外,Huber 等[26]提出了一種分析時間依賴性二氧化碳-鹽水注入的碳封存方法,運用PetraSim 模擬器建立一維流動模型模擬該注入過程,模擬表明,將高密度二氧化碳注入鹽水含水層后進行鹽水驅油,能夠將流動二氧化碳的質量分數降低至10%以下.Lubon 等[27]使用PetraSim 模擬軟件進行了氫氣注入模擬,成功證實當地可以采用合理的氣體回收參數來進行深層含水層的地下儲氫.以上研究表明運用PetraSim 進行數值模擬的可靠性.

本文首先利用圖像法獲取多孔介質中DNAPL污染范圍變化圖像,再利用PetraSim 數值模擬程序獲取DNAPL污染物遷移過程中的飽和度分布圖像.通過將數值模擬得到的圖像與圖像法實驗的結果進行誤差對比分析,評判數值模擬結果準確,通過結合實驗結果與模擬結果分析DNAPL 有機污染物的遷移行為.研究將有助于有效探討飽和多孔介質中DNAPL 的遷移機理與分布擴散特征,為野外三維空間條件下有效的定量評價和預測DNAPL 污染范圍奠定基礎.

1 砂箱實驗

1.1 實驗裝置與材料

室內實驗裝置主要有:二維玻璃砂箱、電子天平、蠕動泵(BT100M)、數碼相機(佳能EOS450D)、相機三角架、燒杯、量筒、溫濕度記錄儀等.

二維砂箱模型由無色鋼化玻璃制成,尺寸為長60cm×寬5cm×高85cm,裝置如圖1所示.在砂箱底部兩側與頂部兩側設置a、b、c、d 四個閥門開關作為進、排水口,用于控制砂箱模型的水位.污染物注入點設置于砂箱頂部中心距頂部15cm 處.填裝背景介質前在砂箱底部鋪設1cm 厚黏土層作為不透水層,保證實驗過程中砂箱整體的密封性[28].

本實驗選用三氯乙烯(TCE)作為DNAPL 典型代表污染物開展實驗,TCE 作為一種無色、易流動的液體,其主要物理化學性質見表1.

為使DNAPL 在二維砂箱中的遷移過程可視化,使用易溶于有機溶液而不溶于水的蘇丹Ⅲ粉末將TCE 溶液染為紅色.利用圖像法拍攝TCE 污染剖面與數值模擬結果進行對比,通過砂箱實驗驗證模擬結果預測DNAPL 遷移行為的可靠性,蘇丹Ⅲ濃度為0.2g/L.實驗選用石英砂作為砂箱內部的多孔介質,石英砂粒徑范圍為40～60 目(0.250～0.425mm),填裝時對石英砂進行粒徑篩選,經過篩選后去除介質內部雜質,并且確保介質的均質性.表2 為物理砂箱實驗中的各項參數.

表2 二維砂箱實驗參數Table 2 Experimental parameters for two-dimensionalsand boxes

1.2 DNAPL 入滲實驗

1.2.1 砂箱裝填 采用40～60 目標準篩網篩取實驗用砂 填裝時進行分層填裝.每次填裝時先使用電子天平稱取石英砂重量,填裝時使用燒杯均勻填入砂礫,填入時以5cm 為一層分層并用木板壓實,使砂層介質分布均勻,均勻填砂至距離砂箱頂部2cm 處,保持砂箱頂部平整并確保砂箱內介質密實.

1.2.2 飽和帶形成 待砂箱裝砂完成后,通過蠕動泵以1mL/s 注水流量將自來水從圖1所示c 號進水閥門注入砂箱內,并將d 號閥門連接一根軟管固定住用以觀察砂箱內部水位變化情況.注水過程中,將砂箱介質孔隙中的空氣完全驅替出來,形成一個飽和二維砂箱模型,待砂箱頂部右側b 閥門有水排出時停止注水.而后將砂箱靜置10h,待水位穩定確保砂箱內處于飽水狀態.

1.2.3 DNAPL 入滲實驗 將染色后的TCE 溶液用注射泵從距砂箱頂部15cm 中心位置定量注入砂箱中,實驗過程中,對污染物注入量、注入流量、排水量進行記錄,當污染物入滲至砂箱底部不透水層時停止注入.注入實驗時長為330min,共注入污染物總量878mL,砂箱排水總量為839mL,存在些許誤差.分析可能原因一是仍然有少量空氣存于裝置中,未及時排出.二是實驗過程中排水量計算為人為估讀,用量筒接取實驗排出水過程中存在殘留水未被估讀在內,造成誤差.

1.2.4 圖像法采集 實驗采用數碼相機記錄實驗過程,實驗開始前,將攝像機固定在砂箱模型前適當位置,確保DNAPL 入滲過程可以被清晰記錄.入滲過程中,每5min 記錄一次,對污染物入滲行為實時拍攝.

1.3 實驗結果與分析

本文通過圖像法獲得DNAPL 在砂箱內的遷移過程圖,如圖2所示.通過繪制DNAPL 在遷移過程中的污染前沿曲線,可視化實驗過程中DNAPL的遷移情況.通過比對各時間節點中DNAPL 入滲圖像,對物理砂箱模型中DNAPL 的遷移規律進行分析.

圖2 砂箱中DNAPL 入滲過程Fig.2 Diagram of DNAPL infiltration process in sandbox

由圖2 可知,污染物的遷移行為主要在重力的影響下進行垂向遷移,同時也在毛細管力的作用下發生橫向遷移.在注入污染物0～60min 時段內,垂向遷移行為最為明顯,垂向遷移長度(24.07cm)遠大于同時段橫向遷移長度(10.92cm),表明該時段內DNAPL 受重力影響更為顯著.90～180min 時段內,DNAPL 以中間較寬,兩端較小的“液滴狀”入滲.水平方向的遷移達到最遠 18.18cm 停止,砂箱中DNAPL 的縱向入滲行為較為連續,并未出現分層現象.在270～330min 時段內DNAPL 呈“指狀”行為入滲,其入滲圖像也表現出不規則特征,DNAPL 的遷移末端也形成不連續殘留,殘留的DNAPL 會占據孔隙體積.

可以看出,垂向上DNAPL 的遷移長度隨污染物的持續注入而逐漸增加,而橫向遷移的范圍在污染物注入180min 后,基本不發生變化.在入滲過程中,DNAPL 橫向遷移的平均速率為0.13cm/min,垂向遷移的平均速率為0.31cm/min,對比可知,污染物垂向遷移的速率明顯大于水平遷移速率.在30min 時,測得砂箱內DNAPL 平面擴散速度為3.80cm2,垂向遷移速度0.59cm/min,橫向遷移速度為0.25cm/min,隨DNAPL 的持續注入,砂箱內DNAPL 的遷移速度逐漸減小,120min 時,測得砂箱中DNAPL 垂向入滲速率為0.35cm/min,橫向遷移速率為0.12cm/min.注入時,DNAPL 垂向速度減小的原因是在DNAPL 的遷移末端發生了殘留現象,隨著遷移時間的變長,殘留的DNAPL 含量越多,孔隙內DNAPL 飽和度增加,其遷移速度也隨之變慢.在330min 時,DNAPL 遷移至砂箱底部,由圖2 可知,DNAPL 在遷移路徑上出現不連續現象,其原因可能在于砂箱的注水飽和過程中介質內部出現優勢通道,在DNAPL 遷移過程中優先通過優勢通道向下入滲,表現為不連續的遷移路徑.實驗過程中DNAPL 污染物總注入量為878mL,在入滲實驗中,飽和多孔介質背景下DNAPL 的污染面積增長率和遷移距離同注入量呈正相關.隨著DNAPL 注入量從30min 時的79.81mL 增加至最后的878mL,垂向遷移距離從最深15.81cm 增加至70cm,橫向遷移距離從最遠 8.12cm 增加至38.86cm,DNAPL 污染面積從 106.88cm2增大至1402.83cm2.飽和多孔介質環境下易于DNAPL 的蓄積,在注入總量不斷增大的情況下,裝置底部DNAPL 不斷蓄積至最大飽和度且圖像法監測的污染面積持續增加.

由實驗結果可以看出,圖像法只能觀測到二維砂箱側面污染物的遷移情況,對于砂箱內DNAPL 入滲過程中污染物的分布情況無法更為直觀的觀測,存在一定局限.為更為精細的研究刻畫出DNAPL 在飽和砂箱中的遷移規律,需建立數值模擬有效刻畫出DNAPL 的遷移與分布規律.

2 數值模擬

2.1 數值模擬方法

本文運用多相流數值模擬軟件PetraSim 建立DNAPL 運移數值模擬模型,通過實驗結果與數值模擬結果的對比,得到最佳的DNAPL 入滲模型參數.本文假定的模型網格按X、Z 方向剖分,模擬區域有效網格設置長60cm,高85cm.網格為不均勻剖分,X方向處左右各設置0.5cm 的隔水邊界,模型底部為零通量邊界,右上角設置排水口,距砂箱頂部中心15cm 處設置污染物泄漏點,概念模型如圖3所示.在建立模型時假定包括:

圖3 數值模型示意Fig.3 Schematic diagram of the conceptual model

(1)模型內基質與流體不可壓縮,水和污染物的密度不變;

(2)不考慮水的蒸發與污染物的揮發;

(3)在DNAPL 模擬入滲前,模型內不存在油相;

(4)忽略降解作用.

相對滲透率函數和毛細管壓力函數是多相流模擬流體在地下水系統中運移分布的重要因子.本文模擬采用三相流相對滲透率函數Stone 模型及毛細管壓力函數Parker 模型[29-30].

相對滲透率函數的具體表達式為:

式中:Kgr、Kwr、 Knr為水相、氣相、NAPL 相的相對滲透率;Sgr、 Swr、 Snr為水相、氣相、NAPL 相殘余飽和度;n 為擬合參數.

毛細管壓力函數具體表達式為:

模型上邊界與大氣相通,壓強為大氣壓強.選取TCE 作為模擬的DNAPL 污染物代表,由于TCE 難溶于水,在水相中溶解度極低,故本文不考慮TCE 的溶解,使其完全以NAPL 相的形式存在模型中.模型溫度為設定為20 ℃.并不考慮熱物理相關參數,模擬時TCE以點源形式從泄漏點開始泄漏.模型的初始條件主要參數見表3.

表3 數值模型中主要參數Table 3 Major parameters in the numerical model

2.2 數值模擬結果與分析

2.2.1 模擬結果 實驗模擬時,DNAPL通過點源的形式從泄漏點開始注入系統中,在模擬注入污染物的前期,DNAPL 以垂向遷移為主,其污染區域飽和度相對較低;隨DNAPL 的持續注入,其污染區域面積也不斷增大,當DNAPL 到達砂槽底部后,在底部聚集并形成長期的飽和態污染源區,并以源區為中心向周圍橫向擴散,飽和度也以源區為核心向外逐漸減小.DNAPL 遷移模擬圖像見圖4.

圖4 DNAPL 飽和度變化Fig.4 Graph of DNAPL saturation variation

從模擬結果可以看出,在模擬至30min 時,污染物在毛細力作用下在泄漏點附近相對均勻的向四周擴散.而在DNAPL 注入60～180min 模擬時段,相對于污染物垂向的遷移,系統兩側基質中DNAPL的橫向入滲明顯滯后,且隨時間推移,系統基質內DNAPL 入滲深度增加,污染范圍增大.在模擬至210min 時,DNAPL 已遷移至模型底部并開始在底部隔水邊界上聚集形成污染池,并且橫向遷移速率開始顯著增大,污染池的長度由19.16cm 增加至270min 時的24.78cm,這表明DNAPL 的運移行徑在毛細力作用下由垂向遷移轉變為以橫向遷移為主.模型在模擬的第330min 時停止DNAPL 的注入,此時在毛細力影響下DNAPL 污染池的長度已從19.16cm 增大至41.00cm,在模型均質介質的背景下,污染池已基本遷移至左右隔水邊界.由圖4 可知,在污染池形成期間,底部源區飽和度隨著DNAPL 的聚集由0.1 增大至0.48,而DNAPL 泄漏點處污染物飽和度降低至0.3,在污染物遷移路徑上的飽和度也隨時間開始降低,其遷移路徑上飽和度范圍由0.38～0.20 降低至0.27～0.12,飽和度隨時間的變化幅度較大.

飽和度模擬結果表明,DNAPL 在含水層中主要受重力作用影響向下遷移,隨時間的增加遷移至底部不透水層,并大量聚集形成DNAPL 污染池.隨著飽和度的持續增加,此時在毛細力的影響下DNAPL 開始以橫向再分布遷移為主,橫向遷移距離與污染分布面積持續增大,最終在模型飽水帶中達到油-水平衡狀態.在 DNAPL 遷移至含水層底部過程中,部分DNAPL 受介質顆粒的吸附作用長期殘留在顆粒孔隙中,與污染池一起對地下水形成長期污染.

2.2.2 數據分析 為將數值模擬所獲結果與實驗實際結果進行對比,在實驗結果與模擬結果的二維平面圖像Z=0.6m 與X=0.3m 處分別做實驗結果橫剖面線A-A1、模擬結果橫剖面線A-A2 和實驗結果縱剖面線B-B1 與模擬結果縱剖面線B-B2 進行污染物遷移距離對比分析,綜合物理模型與數值模型不同時刻DNAPL 的遷移入滲情況分析DNAPL 在含水層中的運移規律.所選取的剖面線示意圖如圖5所示.二維砂箱實驗與模擬中DNAPL 的相關運移參數與各時刻的相對誤差見表4.相對誤差計算公式為:

表4 DNAPL 運移特征相關參數Table 4 Parameters related to DNAPL transport characteristics

圖5 剖面線示意Fig.5 Schematic diagram of the profile line

式中:x 為模擬值,cm;μ為實驗值,cm.

由表4 可以看出,模擬與實驗垂向上DNAPL的遷移長度隨污染物的持續注入而逐漸增加,但在210min 后,DNAPL 遷移至底部隔水層,其垂向污染范圍基本不再變化,且無論是模擬實驗或是物理實驗,在污染物注入過程中,DNAPL 垂向遷移距離始終大于水平方向的遷移距離.室內砂箱實驗與數值模擬得到的DNAPL 入滲過程對比圖如圖6所示,在0～30min 注入DNAPL 的初期,DNAPL 在注入點附近的擴散主要受毛細力的影響,其垂向的遷移距離與橫向遷移距離相對平均.在 30min 之后,DNAPL 運移行為主要受重力影響,垂向遷移的距離遠大于橫向距離,觀察圖6 模擬時的飽和度分布變化圖可知,以注入點為中心,周圍基質內油飽和度開始不斷增加,其分布變化趨勢與圖像法的分析結果相一致.在實驗與模擬時間達到 210min 時,DNAPL 遷移至底部不透水層,且在重力作用下,DNAPL 由于垂向方向的遷移,污染區域發生明顯的下移,分析是由于DNAPL 開始在不透水層表面聚集形成污染池,DNAPL 飽和度增大.隨著污染池的形成,DNAPL 的遷移在毛細力作用下以橫向再分布遷移為主.分析210～330min 時的對比圖時,DNAPL 向左右兩側發生橫向遷移,其底部飽和度也不斷增大,而遷移路徑上的污染物在重力作用下運移至底部污染池,可看出油飽和度明顯降低,但由于介質顆粒的吸附作用,仍有少量污染物殘留在孔隙空間內.再分布期間污染物分布的模擬結果與實驗結果基本一致.

圖6 均質砂箱中實驗圖像與模擬結果對比Fig.6 Comparison of experimental images and simulation results in a homogeneous sand box

為進一步評估數值模擬方法可視化DNAPL 源區遷移過程的有效性,通過相對誤差分析各時刻模擬值與實驗值的橫向遷移距離、垂向遷移距離與污染面積的誤差.通過繪制實驗值與模擬值的關系曲線作為整組數據分析的誤差指標,采用擬合決定系數R2來體現模擬值與實驗值整體間的差異,R2的取值在[0,1]之間,越接近1,說明回歸擬合效果越好,模擬值與實驗值整體差異越小.

擬合結果表征了模擬結果與實驗結果中污染物垂向分布范圍、橫向分布范圍與二維平面污染物面積的差異,表征結果如圖7所示.模擬結果的二維污染物平面擴散趨勢與實驗結果非常接近.由圖7可知,模擬結果與實驗結果中污染物垂向范圍的相對誤差范圍在-2.96%～5.95%,而橫向遷移前90min的相對誤差相對于垂向遷移的更大,最大可達23.41%,這可能是由于橫向遷移行為受毛細管力影響大,實驗時泄漏點注入污染物時附近的孔隙空間有所改變,而隨著注入時間的增加,實驗系統達到穩定,相對誤差逐漸降低.污染面積分布的相對誤差范圍在0.03%～19.39%,其中30、90 和270min 的相對誤差較大,高于15%,這是由于數值模擬存在的局限性,無法如圖像法那般精細的刻畫出二維砂箱中污染物的展布范圍.表征污染物各時刻橫向范圍、垂向范圍與污染面積模擬值與實驗值差異性的擬合決定系數R2分別為0.8827、0.9919、0.9832.可以看出,模擬結果展現出與實際情況極大的相似性,結果表明所建立的數值模型能準確模擬出室內實驗結果.

圖7 DNAPL 入滲過程橫向、垂向遷移距離和污染面積關系曲線圖與各時刻相對誤差Fig.7 Horizontal、vertical migration distance and contaminated area during DNAPL infiltration process and relative error diagram at each moment

3 結論

3.1 DNAPL 污染物在飽和均質砂土介質遷移過程中主要受重力和毛細力作用,污染物在入滲前期主要受重力影響呈“液滴”形垂向擴散.在垂向擴散至底部不透水層后,DNAPL 開始聚集形成污染池,此時DNAPL 遷移行為轉變為由毛細力驅動引導的橫向再分布擴散,最后在底部形成長期的“池狀”污染源區.

3.2 基于Petrasim 程序模擬驗證飽和多孔介質物理模型實驗的DNAPL 入滲過程,探討DNAPL 源區運移分布的形成過程.分析圖像法與數值模擬形成的污染物擴散面積誤差,橫向擴散范圍、垂向擴散范圍與污染擴散面積的相對誤差范圍分別為0.04%～23.41%、-2.96%～5.95%和0.03%～19.39%.

模擬值與實驗值的擬合決定系數R2分別為0.8827、0.9919、0.9832.模擬結果與實驗值接近,模擬效果良好,說明Petrasim 程序探究飽和多孔介質中DNAPL 運移行為與污染展布范圍是可行性的,在DNAPL 污染范圍的定量評價和預測方面具有一定的應用價值.